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Part  II  of  this  thesis  examines  query-retrieval  problems  concerning  distributions 
of  points  in  euclidean  space.  In  this  part,  we  describe  a  new  technique  for  solving  a 
variety  of  query-retrieval  problems  in  optimal  time  with  optimal  or  near-optimal  space 
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for  circular  range  searching,  half-space  range  searching,  and  computing  b-nearest 
neighbors  in  a  variety  of  metrics.  For  each  problem  and  each  query,  the  response 
to  the  query  is  provided  in  O(k)  or  0(k  +  logn)  time  where  k  is  the  size  of  the 
response  and  n  is  the  size  of  the  problem.  (E.g.,  for  the  n-point  k- nearest  neighbors 
problem,  the  nearest  neighbors  of  any  query  point  are  provided  in  0(k  +  logn) 
steps.)  Depending  on  the  problem  being  solved,  the  space  required  for  the  data 
structure  is  either  linear  or  O(nlogn).  Hence,  the  time  bounds  are  optimal  and 
the  space  bounds  are  optimal  or  near-optimal.  Previously  known  data  structures 
for  these  problems  required  a  factor  of  n(logn(loglogn)2)  or  fl(logn  log  logn)  more 
space  and/or  more  time  to  answer  each  query. 

Our  compaction  technique  incorporates  planar  separators,  filtering  search,  and 
the  probabilistic  method  for  discrepancy  problems.  The  fundamental  idea  is  that  k,K- 
order  Voronoi  diagrams  (and  other  suitable  proximity  diagrams)  can  be  compacted 
from  kon]n  space  to  0(n )  space  and  still  retain  all  the  information  that  is  essential 
for  solving  query  problems. 


Results  in  Computational  Geometry: 
Geometric  Embeddings 

and  Query-Retrieval  Problems 

bv 

Mark  David  Hansen 

A.B.,  Mathematics 
Cornell  University 
(1986) 

M.S..  Mathematics 
University  of  Chicago 
(1987) 

Submitted  to  the  Department  of  Mathematics 
in  partial  fulfillment  of  the  requirements  for  the  degree  of 

Doctor  of  Philosophy 

at  the 

MASSACHUSETTS  INSTITUTE  OF  TECHNOLOGY 

January  1991 

©  Massachusetts  Institute  of  Technology  1991 


Signature  of  Author 


Certified  by 


Accepted  by 


Department  of  Mathematics 
January  11,  1991 


F.  Thomson  Leighton 
Professor  of  Applied  Mathematics 
Thesis  Supervisor 


Daniel  J.  Kleitman.  Chairman 
Applied  Mathematics  Committee 


Accepted  by 


Sigurdur  Helgason,  Chairman 
Departmental  Graduate  Committee 


DTI« 


Results  in  Computational  Geometry: 


Geometric  Embeddings 
and  Query-Retrieval  Problems 

by 

Mark  David  Hansen 


Submitted  to  the  Department  of  Mathematics 
on  January  11,  1991.  in  partial  fulfillment  of  the 
requirements  for  the  degree  of 
Doctor  of  Philosophy 


Abstract 


Accccion  For  j  j 

NT  IS  CRA&I 

C  i  iC  TAB 

n 

U.  j  uounced 

□ 

Justification 


By 


Oi-i  id-. 

.tijn/ 

Availability  Codes 

' 

r,  V 

AvSil  «. 
Spe 

md  i  or 
ciol 

Many  fundamental  questions  in  computational  geometry  arise  from  the  consider¬ 
ation  of  distributions  of  points  in  euclidean  space.  This  thesis  explores  two  important 
areas  of  computational  geometry  in  this  setting:  geometric  embeddings  and  query- 
retrieval  problems.  Each  area  is  addressed  in  a  separate  part  of  the  thesis. 

Part  I  examines  the  geometric  embedding  problem  for  many  of  the  graphs  which 

are  important  in  the  study  of  parallel  computation ,[33].  Given  an  undirected  graph 

( _  •-  . . .  " ' 

G  with  n  vertices,  and  a  set  P  of  n  points  in  the  plane,  the  geometric  embedding 

problem  consists  of  finding  a  bijection  from  the  vertices  of  G  to  the  points  in  the  plane 
which  minimizes  the  sum  total  of  edge  lengths  of  the  embedded  graph.  We  give  fast 
approximation  algorithms  for  embedding  d-dimensional  grids  in  the  plane  which  are 
within  a  factor  of  O(logn)  times  optimal  cost  for  d  >  2  and  0(log2  n)  for  d  =  2.  We 
also  show  that  any  embedding  of  a  hypercube,  butterfly,  or  shuffle-exchange  graph 
must  be  within  an  <3(log  n)  factor  of  optimal  cost.  When  the  points  of  P  are  randomly 
distributed,  or  arranged  in  a  grid,  we  give  a  polynomial  time  algorithm  which  can 
embed  arbitrary  weighted  graphs  in  these  points  with  cost  within  an  0( log2  n)  factor 
of  optimal.  Many  of  these  results  extend  to  higher  dimensions  when  P  C  Rd- 
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Aside  from  the  intrinsic  mathematical  interest  of  these  problems,  they  also  have 
applications  in  the  field  of  parallel  processing.  For  example,  we  show  how  the  al¬ 
gorithms  which  we  develop  for  geometric  embeddings  can  be  used  to  give  solutions 
which  are  within  an  0(log 2  yV)  factor  of  optimal  to  problems  of  performance  opti¬ 
mization  fjr  array-based  parallel  processors  in  the  following  areas:  communication 
load  balancing,  dynamic  allocation  of  jobs  to  processors,  reconfiguring  around  faults, 
and  simulating  other  architectures. 

'  Part  II  of  this  thesis  examines  query-retrieval  problems  concerning  distributions 
of  points  in  euclidean  space.  In  this  part,  we  describe  a  new  technique  for  solving  a 
variety  of  query- retrieval  problems  in  optimal  time  with  optimal  or  near-optimal  space  * 
[2].  In  particular,  we  use  the  technique  to  construct  algorithms  and  data  structures 
for  circular  range  searching,  half-space  range  searching,  and  computing  fc-nearest 
neighbors  in  a  variety  of  metrics.  For  each  problem  and  each  query,  the  response 
to  the  query  is  provided  in  O(k)  or  0(k  -f  logn)  time  where  k  is  the  size  of  the 
response  and  n  is  the  size  of  the  problem.  (E.g.,  for  the  n-point  fc-nearest  neighbors 
problem,  the  k- nearest  neighbors  of  any  query  point  are  provided  in  0(k  -f  logo) 
steps.)  Depending  on  the  problem  being  solved,  the  space  required  for  the  data 
structure  is  either  linear  or  O(nlogn).  Hence,  the  time  bounds  are  optimal  and 
the  space  bounds  are  optimal  or  near-optimal.  Previously  known  data  structures 
for  these  problems  required  a  factor  of  f!(log  n(log  log  n)2)  or  D(log  n  log  log  n)  more 
space  and/or  more  time  to  answer  each  query. 

Our  compaction  technique  incorporates  planar  separators,  filtering  search,  and 
the  probabilistic  method  for  discrepancy  problems.  The  fundamental  idea  is  that  kth- 
order  Voronoi  diagrams  (and  other  suitable  proximity  diagrams)  can  be  compacted 
from  k°^'n  space  to  0(n)  space  and  still  retain  all  the  information  that  is  essential 
for  solving  query  problems. 

Thesis  Supervisor:  F.  Thomson  Leighton 
Title:  Professor  of  Applied  Mathematics 
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Introduction 


This  thesis  examines  problems  in  computational  geometry  which  arise  from  the  con¬ 
sideration  of  distributions  of  points  in  euclidean  space.  The  fundamental  case  concerns 
a  set  of  n  points  P  =  px ,p2,  •  •  •  ,pn  in  the  euclidean  plane.  The  class  of  problems  we 
present  involves  constructing  the  optimal  embedding  of  of  a  given  graph,  such  as  a 
mesh,  in  P.  Such  problems  arise,  for  example,  when  one  tries  to  determine  how  to 
optimally  link  together  a  group  of  distributed  computers  into  a  mesh  or  hypercube 
network.  This  class  of  problems  is  referred  to  as  the  geometric  embedding  problem. 
In  general,  geometric  embedding  problems  are  NP-complete.  In  this  thesis,  we  give 
approximation  algorithms  which  yield  nearly  optimal  solutions  to  a  large  variety  of 
geometric  embedding  problems  in  polynomial  time. 


The  second  class  of  geometry  problems  we  consider  are  query-retrieval  problems. 
In  this  case,  we  are  interested  in  “retrieving”  the  answers  to  a  set  of  queries  about 
P.  For  example:  “Among  the  n  points  in  the  plane,  find  the  k  points  that  are  closest 
to  position  (x,,p,)  for  i  =  1  . . .  m.”  The  task  is  to  preprocess  P  and  devise  a  data 
structure  so  that  each  query  can  be  answered  as  quickly  as  possible.  In  this  thesis,  we 
are  primarily  concerned  with  minimizing  the  space  used  by  the  query- retrieval  data 
structure,  while  maintaining  optimal  time  answers  to  individual  queries. 
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INTRODUCTION 


Geometric  Embeddings 

Part  I  of  this  thesis  describes  the  geometric  embedding  results.  In  Chapter  1  a  precise 
definition  of  the  geometric  embedding  problem  is  given,  the  major  contributions  of 
this  thesis  in  that  area  are  described,  and  applications  of  this  work  to  solving  par¬ 
allel  processing  problems  are  discussed.  Chapters  2  and  3  describe  approximation 
algorithms  for  near  optimal  embeddings  of  graphs  important  to  the  theory  of  parallel 
computation:  the  hypercube,  butterfly,  shuffle-exchange  graph,  and  mesh.  Chapter  4 
presents  another  class  of  approximation  algorithms  for  geometric  embeddings  which 
are  based  on  a  powerful  separator  theorem  due  to  Leighton  and  Rao  [38].  These 
algorithms  give  near-optimal  embeddings  of  any  graph  in  the  plane,  provided  some 
assumptions  are  made  about  the  distribution  of  points  P.  Finally,  in  Chapter  5  we 
describe  in  detail  the  relationship  between  the  embedding  results  and  problems  in 
parallel  computation. 


Query- Retrieval  Problems 

Part  II  of  this  thesis  describes  the  query- retrieval  results.  In  Chapter  6  we  give  some 
background  on  the  importance  of  data  structures  in  the  solution  of  query-retrieval 
problems.  The  main  applications  to  ^-nearest  neighbor  search,  circular  range  search, 
and  half-space  range  search  are  then  described.  Chapter  7  contains  a  description  of 
the  fundamental  compaction  technique  which  allows  us  to  greatly  reduce  the  space 
required  to  store  a  Voronoi  diagram’s  essential  information  for  solving  query- retrieval 
problems.  Lastly,  in  Chapter  8,  the  application  of  this  compaction  technique  to  the 
three  problems  described  in  Chapter  6  is  presented  in  detail. 


Part  I 

Geometric  Embeddings 
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Chapter  1 


Overview  of  Geometric 
Embedding  Results 


1.1  Geometric  Embeddings 

A  geometric  embedding  of  G  =  (V,  E)  into  a  set  of  points  P  in  euclidean  space  is 
a  bijection  /  :  V  — ►  P.  We  are  interested  in  finding  /  which  minimizes  the  total 
edge  length  of  the  graph  induced  on  P:  J2(u,v)eE  dist(f(u),  f(v)).  Here,  dist  is  the 
euclidean  metric.  The  most  famous  example  of  a  geometric  embedding  problem  is 
the  euclidean  traveling  salesman  problem.  In  this  case,  G  is  a  cycle  on  N  nodes  and 
many  approximation  algorithms  for  finding  the  near-optimal  cost  embedding  of  G 
into  points  in  R 2  are  known.  One  example  of  such  an  approximation  algorithm  is 
due  to  Christofides  [23].  In  this  paper,  we  will  consider  the  geometric  embedding 
problem  for  graphs  which  are  important  in  the  theory  of  parallel  computation:  grids, 
hypercubes,  butterflies,  and  shuffle  exchange  graphs. 

We  will  also  consider  some  special  cases  of  the  geometric  embedding  problem  in 
the  plane  and  higher  dimensional  spaces  when  the  set  of  points  P  are  distributed 
randomly,  or  arranged  in  a  grid.  With  these  restrictions  on  P  we  are  able  to  address 
the  more  general  problem  of  embedding  weighted,  undirected  graphs  in  Rd.  For  a 
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weighted  graph,  our  goal  is  to  minimize  the  weighted  sum  of  edge  lengths  in  the 
graph  induced  on  P:  w{u,v)dist(f(u),  f(v)).  Here,  w(u,v)  is  the  weight  of 

edge  (u.u)  €  E.  In  this  case  we  may  also  consider  dist  to  be  the  manhattan  (L i) 
metric.  Some  recent  results  in  multicommodity  flow  [38]  allow  us  to  achieve  near- 
optimal  embeddings  of  arbitrary  weighted  graphs  in  these  two  special  cases.  All  of 
these  results  have  important  applications  to  parallel  processing. 


1.2  Mathematical  Results 

Polynomial  time  algorithms  for  finding  optimal  or  0(logiV)  times  optimal  embed¬ 
dings  in  Rd  of  simple  structures  such  as  cycles,  trees,  or  stars  are  familiar  to  many 
researchers  [11,  23].  In  this  paper,  we  study  the  problem  of  embedding  more  com¬ 
plicated  graphs  in  euclidean  space.  The  structures  we  study  include  many  of  the 
important  graphs  from  the  theory  of  parallel  computation:  hypercubes,  butterflies, 
shuffle-exchange  graphs,  meshes,  cubes,  and  higher  dimensional  grids. 

A  straightforward  application  of  known  routing  results  proves  that  any  embedding 
of  the  high-bandwidth,  low  diameter  hypercube-like  graphs  is  within  an  0(log  N)  fac¬ 
tor  of  optimal.  Meshes,  cubes,  and  higher-dimensional  grids  do  not  have  the  nice  rout¬ 
ing  properties  of  these  graphs  and  therefore  embedding  them  efficiently  constitutes  a 
significantly  more  difficult  problem.  Section  3  presents  the  analytic  techniques  which 
enable  us  to  give  a  fast  algorithm  for  embedding  these  graphs  within  an  0(log  N)  or 
0(log2  N)  factor  of  optimal.  Some  of  these  techniques  can  be  generalized  and  used 
together  with  separator  approximation  algorithms  [38]  to  give  near-optimal  embed¬ 
dings  of  arbitrary  weighted  graphs  into  uniform  distributions  of  points  in  euclidean 
space. 

The  embedding  results  in  this  paper  are  grouped  into  two  categories  based  on  the 
arrangement  of  the  set  of  points  P  into  which  we  are  embedding  the  graph  G.  P 
could  be  (1)  an  arbitrary  set  of  points,  or  (2)  a  array  of  points  or  a  random  set  of 
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uniformly  distributed  points. 

1.2.1  Arbitrary  Points 

For  hypercubes,  butterflies,  and  shuffle-exchange  graphs  it  is  straightforward  to  show 
that  any  embedding  is  near-optimal.  In  fact,  in  this  case,  we  need  not  even  assume 
that  the  points  P  are  in  the  plane,  but  only  that  they  are  in  some  graph  which  obeys 
the  triangle  inequality. 

Theorem  1.2.1  Given  a  graph  G  which  is  an  N -node  hypercube,  butterfly,  or  shuffle- 
exchanqe  graph,  and  a  set  P  of  points  connected  by  edges  which  obey  the  triangle 
inequality,  any  embedding  of  G  into  P  has  cost  within  an  0(\ogN)  factor  of  optimal. 


For  grids,  cubes,  and  higher  d-dimensional  grids,  the  problem  is  substantially  more 
difficult,  and  we  give  a  fast  algorithm  for  near  optimal  embeddings  in  the  plane. 

Theorem  1.2.2  There  exists  an  algorithm  which,  when  given  a  d-dimensional  grid  G 
on  N  nodes  and  a  set  P  of  N  points  in  the  plane,  in  0(N  log  N)  time  finds  a  geometric 
embedding  with  cost  a  factor  of  0( log2  N)  times  optimal  when  G  is  a  square  mesh 
(d  =  2 )  and  0(\ogN)  times  optimal  when  G  is  a  cube  or  higher  dimensional  grid 

(d>  2). 

In  fact,  Theorem  1.2.2  generalizes  to  the  case  P  C  Rk  in  the  following  manner. 
We  can  embed  grids  of  dimension  k  into  Rk  within  an  0(log2  N)  factor  of  optimal 
cost,  and  grids  of  dimension  d  >  k  within  an  0(log  N)  factor. 

1.2.2  Arrays  of  Points  and  Randomly  Distributed  Sets  of 
Points 

We  can  achieve  significantly  more  general  embedding  results  when  we  place  some 
restrictions  on  the  arrangement  of  the  points  P.  If  P  is  a  set  of  evenly  spaced  points 
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arranged  in  a  d-dimensional  array  in  Rd  for  some  fixed  constant  d,  then  we  can  use 
the  graph  separator  approximation  results  of  Leighton  and  Rao  [38]  to  construct  a 
polynomial  time  algorithm  which  produces  a  geometric  embedding  of  an  arbitrary 
weighted  graph  G  into  P  with  cost  within  an  0(log2  N)  factor  of  optimal. 

Theorem  1.2.3  There  exists  a  polynomial  time  algorithm  which,  when  given  an  ar¬ 
bitrary  N -node  weighted  graph  G  and  a  d-dimensional  array  of  N  points  P  in  Rd. 
embeds  G  in  P  with  cost  within  an  <9(log2  N)  factor  of  optimal. 

The  same  theorem  holds,  with  high  probability,  when  P  is  a  random  set  of  uni¬ 
formly  distributed  points  in  Rd. 

Theorem  1.2.4  There  exists  a  polynomial  time  algorithm  which,  when  given  an  ar¬ 
bitrary  N -node  graph  G  and  a  random  set  of  N  uniformly  distributed  points  P  in  Rd, 
with  high  probability  embeds  G  in  P  with  cost  within  an  £>(log2  N)  factor  of  optimal. 


1.3  Applications 

1.3.1  Efficient  Parallel  Computation 

Efficient  parallel  processing  requires  effective  algorithms  for  embedding  large  compu¬ 
tational  problems  into  parallel  architectures.  Given  a  computation  graph  for  a  parallel 
algorithm  which  consists  of  process  nodes  and  edges  between  nodes  whose  processes 
exchange  data,  and  a  fixed  parallel  machine  architecture,  we  would  like  to  embed  this 
computation  graph  into  the  architecture  in  a  manner  which  makes  efficient  use  of  the 
machine’s  interconnect  network  for  inter-process  communication.  A  poor  embedding 
could  result  in  an  assignment  of  process  nodes  to  processors  which  clogs  the  network 
with  inter-process  communication  and  greatly  degrades  running  time. 

We  measure  the  efficiency  of  an  embedding  of  a  computation  graph  into  a  parallel 
architecture  in  terms  of  the  communication  load  this  embedding  induces  on  the  ma- 
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chine.  Communication  load  is  a  measure  of  the  total  volume  of  traffic  an  algorithm 
produces  on  a  network,  and  is  defined  to  be  the  sum  total  of  distances  that  inter¬ 
process  messages  travel  in  the  interconnect  network.  As  a  corollary  of  Theorem  1.2.3 
we  show  how  to  embed  any  computation  graph  in  an  array-based  parallel  machine 
with  communication  load  within  an  0(log2  N )  factor  of  optimal.  This  corollary  can 
also  be  interpreted  to  mean  that  we  can  reconfigure  an  array-based  machine  to  sim¬ 
ulate  any  other  architecture  within  an  0(log2  N)  factor  of  optimal  communication 
load. 

Theorems  1.2.1,  1.2.2,  and  1.2.4  allow  us  to  give  algorithms  for  the  efficient  dy¬ 
namic  allocation  of  resources  on  a  multiprocessor.  Suppose  that  our  machine  has  a 
number  of  users  who  are  continually  submitting  jobs  for  parallel  processing.  The  ma¬ 
chine’s  operating  system  dynamically  allocates  processors  to  jobs,  and  as  jobs  finish. 
holes  of  idle  processors  open  up  in  the  network.  For  arbitrary  holes  in  an  array-based 
processor,  our  results  give  algorithms  for  embedding  jobs  with  near-optimal  commu¬ 
nication  load  when  these  jobs  have  common  computation  graphs  such  as  grids,  cubes, 
hypercubes,  shuffle-exchange  graphs,  or  butterflies.  When  the  holes  are  assumed  to 
be  randomly  and  uniformly  distributed,  we  can  give  efficient  embedding  algorithms 
for  any  computation  graph. 

1.3.2  Network  Reconfigurat’on  and  Wafer-Scale  Integra¬ 
tion 

Consider  a  nation-wide  network,  similar  to  one  operated  by  the  weather  service,  of 
computation  centers  linked  over  standard  telecommunication  channels.  If  the  network 
is  tracking  a  storm  system,  it  will  need  to  be  continually  reconfigured  as  centers  enter 
and  drop  out.  Theorems  1.2.1,  1.2.2,  and  1.2.4  provide  algorithms  for  dynamically 
reconfiguring  in  a  manner  which  produces  networks  with  near-optimal  communication 
overhead. 

In  wafer-scale  integration  [37]  reconfiguring  the  live  cells  on  a  silicon  wafer  into 
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a  usable  network  comprises  a  geometric  embedding  problem.  Similarly,  consider  re¬ 
configuring  the  live  processors  on  an  array-based  machine  with  faults.  Our  results 
indicate  how  to  reconfigure  for  near-optimal  communication  load  on  the  water  or  ar¬ 
ray  in  the  presence  of  arbitrary  faults  when  the  embedded  network  is  a  grid,  cube,  or 
hypercube-like  graph.  If  faults  are  random  or  uniformly  distributed,  we  can  reconfig¬ 
ure  well  for  any  embedded  network. 


1.4  Related  Work 

Several  special  cases  of  finding  optimal  solutions  to  the  geometric  embedding  problem 
have  been  proven  A'P-hard.  These  include  the  cases  when  G  is  a  traveling  salesman 
tour  or  binary  tree  [32].  Christofides!  approximation  algorithm  for  the  geometric 
traveling  salesman  achieves  embeddings  which  are  1.5  times  optimal.  Papadimitriou 
and  Yanakakis  [40]  have  further  shown  that  the  geometric  embedding  problem  is 
N P-hard  for  large  classes  of  trees. 

In  a  recent  paper,  Bern,  Karloff,  Raghavan,  and  Schieber  [11]  studied  geometric 
embeddings  of  binary  trees,  cycles,  and  stars.  They  present  a  number  of  interesting 
results  including  very  fast  approximation  algorithms  for  embedding  these  graphs  in 
the  line  and  plane.  The  focus  of  their  paper  is  on  finding  fast  algorithms  for  efficiently 
embedding  these  simple  structures.  In  contrast,  the  present  paper  concentrates  on 
finding  near-optimal  embedding  techniques  for  some  of  the  more  complex  graphs 
studied  in  the  theory  of  parallel  processing.  In  the  case  of  d-dimensional  grids,  the 
approximation  algorithm  we  present  also  happens  to  be  very  fast. 

1.5  Outline  of  Part  I 

The  remainder  of  Part  I  is  divided  into  four  chapters.  Chapter  2  presents  Theorem 
1.2.1  as  a  straightforward  corollary  of  some  well-known  routing  results.  Here  we 
indicate  the  connection  between  routing  results  and  geometric  embeddings  which 
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will  be  used  in  the  analysis  of  the  grid  embedding  algorithm.  Chapter  3  presents 
the  approximation  algorithm  and  mathematical  analysis  necessary  to  prove  Theorem 
1.2.2.  Chapter  4  contains  the  probabilistic  analysis  and  approximation  algorithms 
for  Theorems  1.2.3  and  1.2.4.  Finally.  Chapter  5  discusses  some  applications  of  these 
theorems  to  parallel  processing. 
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Chapter  2 


Hypercube  Related  Graphs 


In  this  chapter  we  use  well  known  routing  results  to  prove  Theorem  1.2.1  and  show 
that  any  embedding  of  the  high-bandwidth,  low-diameter  hypercube-like  graphs  is 
within  an  0(log  N)  factor  of  optimal.  Notice  that  the  proof  given  below  does  not  use 
any  properties  of  euclidean  space,  but  only  requires  that  the  metric  being  embedded 
into  obey  the  triangle  equality. 

Theorem  1.2.1  Given  a  graph  G  which  is  an  Ar-node  hypercube,  butterfly,  or 
shuffle-exchange  graph,  and  a  set  P  of  points  connected  by  edges  which  obey  the 
triangle  inequality,  any  embedding  of  G  into  P  has  cost  within  an  O(logiV)  factor  of 
optimal. 

Proof:  Let  P  be  a  complete  graph  on  N  points  with  edge  weights  which  obey 
the  triangle  inequality.  (P  being  points  in  the  plane  with  the  euclidean  metric  is 
hence  a  special  case.)  For  any  subset  E  of  edges  in  P ,  let  ||£’||  be  the  sum  of  the 
weights  of  these  edges. 

First  consider  the  jV-node  hypercube  H.  Let  fopt  :  H  — ►  P  be  an  optimal  em¬ 
bedding  of  a  hypercube  in  P.  Then  let  g  :  H  — *  P  be  any  other  embedding.  Define 
f0Vt{H)  and  g{H)  to  be  the  hypercube  subgraphs  of  P  induced  by  the  edges  of  H. 
Group  the  edges  of  H  by  dimension,  so  that  we  have  groups  E\,  ...,  £|gJv  where  the 
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edges  in  Et  connect  points  across  the  ith  dimension.  Then  g(E,)  defines  a  matching  in 
P.  The  following  well-known  lemma  due  to  Benes  [9]  indicates  that  we  can  connect 
up  the  pairs  in  this  matching  via  paths  in  fopt{H).  using  each  edge  of  fopt(H)  at  most 
0(1)  times: 

Lemma  2.0.1  (Benes)  Any  permutation  of  A’  nodes  in  the  N-nnde  hypercube  H 
can  be  routed  so  that  each  edge  of  H  is  used  at  most  0(1)  times. 

Hence,  by  the  triang'e  inequality,  ||<;(.E,)||  <  0(||/opt(i/)||).  Summing  over  the 
logiV  dimensions,  we  see  that  j|^(/f)||  <  0(lg  iV)||/opt(/f)||. 

Now  let  B  be  the  A'-node  butterfly.  Again,  let  fopt  :  B  — ►  P  be  an  optimal 
embedding  of  a  butterfly  in  P  and  let  g  :  B  — ►  P  be  any  other  embedding.  Then  the 
edges  of  g(B)  can  be  decomposed  into  4  disjoint  partial  matchings.  (B  has  degree 
4).  Another  familiar  routing  lemma  [9]  shows  that  each  matching  can  be  routed  in 
fopt(B )  using  edges  at  most  O(logiV)  times: 

Lemma  2.0.2  (Benes)  Any  permutation  of  N  nodes  in  the  N -node  butterfly  B,  or 
N -node  shuffle  exchange  graph  S,  can  be  routed  so  that  each  edge  of  B  or  S  is  used 
at  most  0( log  N)  times. 

Hence,  using  the  triangle  inequality  ||^(B)||  =  0(\og  N)\\fopt(B)\\.  The  proof  for 
the  shuffle  exchange  graph  is  similar.  □ 

Grids  do  not  have  the  high  bandwidth  and  low  diameter  of  the  three  graphs  studied 
above.  Hence  we  cannot  expect  to  get  an  O(logyV)  factor  approximation  algorithm 
for  embedding  grids  using  these  simple  routing  ideas  alone.  We  will,  however,  use 
routing  results  in  the  analysis  of  our  approximation  algorithm  for  grids. 


Chapter  3 

d-Dimensional  Grids 


3.1  Overview 


Chapter  3  is  devoted  to  the  proof  of  Theorem  1.2.2. 


Theorem  1.2.2  There  exists  an  algorithm  which,  when  given  a  d-dimensional  grid  G 
on  N  nodes  and  a  set  P  of  N  points  in  the  plane,  in  0(N  log  N)  time  finds  a  geometric 
embedding  with  cost  a  factor  of  0(log2  N)  times  optimal  when  G  is  a  square  mesh 
( d  =  2)  and  0(\ogN )  times  optimal  when  G  is  a  cube  or  higher  dimensional  grid 

(d>  2). 

We  begin  by  using  some  routing  results  on  grids  to  prove  that  certain  sub-regions 
of  the  plane  are  always  densely  packed  with  points.  A  divide  and  conquer  subroutine 
is  then  invoked  to  embed  pieces  of  the  grid  into  these  dense  regions.  This  subroutine 
is  proven  to  give  nearly  optimal  embeddings  of  these  pieces  provided  the  image  points 
are  contained  in  a  sufficiently  dense  region.  Finally,  we  show  that  these  pieces  can  be 
hooked  together  to  form  an  embedding  of  the  entire  grid  with  near-optimal  cost. 
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3.2  A  Review  of  Routing  Results  for  Grids 

Within  constant  factors,  the  best  routing  result  which  can  be  achieved  for  d- 
dimensional  grids  states  that: 

Lemma  3.2.1  Any  permutation  of  an  m  point  subset  of  a  d-dimensional  gnd  G  can 
be  routed  using  each  edge  of  G  at  most  2|’ma'j  times. 

In  order  to  prove  Lemma  3.2.1,  we  begin  by  proving  a  somewhat  less  general  result 
concerning  the  congestion  required  to  route  permutations  of  an  entire  d-dimensional 
grid: 

Lemma  3.2.2  Any  permutation  of  the  N  =  nd  points  of  a  d  dimensional  grid  G  can 
be  routed  using  each  edge  at  most  2 N?  times. 

Proof:  The  proof  proceeds  by  induction  on  d.  For  d  =  1  we  have  a  line,  and  any 
greedy  routing  will  work.  Now  suppose  that  the  theorem  holds  for  d  —  1.  Pick  a 
dimension  and  consider  G  to  be  composed  of  n  d  —  1  dimensional  planes  connected 
by  columns  along  that  dimension.  Each  point  is  contained  in  a  unique  column  with 
n  —  1  other  points.  Consider  each  source  and  destination  pair  of  the  permutation. 
By  induction  we  will  be  finished  if  we  can  permute  points  within  their  columns  so 
that  each  source  and  destination  pair  end  up  in  the  same  plane.  For  at  this  stage 
we  would  then  have  n  similar  routing  problems  on  the  n  planes.  Each  edge  of  the 
columns  in  our  chosen  dimension  would  have  been  used  at  most  2n  =  2 N?  times,  and 
later  stages  in  the  routing  never  use  that  dimension  again.  So  it  suffices  to  show  that 
we  can  permute  the  points  within  a  column  so  source  and  destination  points  always 
end  up  in  the  same  plane. 

Consider  a  bipartite  graph  with  one  node  on  each  side  for  each  column.  Draw  an 
edge  for  each  source  and  destination  pair  between  any  two  columns.  See  Figure  3-1. 

The  result  is  a  bipartite  graph  on  2n  nodes  which  is  n-regular.  It  is  well  known 
that  such  a  graph  can  be  n-colored.  That  is,  color  the  edges  of  the  graph  with  n 
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Figure  3-1:  Columns  connected  by  source  and  destination  pairs  yield  an  n-regular 
graph  which  can  be  n-colored. 

colors  so  that  each  node  has  no  two  edges  leaving  it  of  the  same  color.  Since  each 
colored  edge  corresponds  to  a  source  and  destination  pair,  we  can  assign  each  pair  to 
the  plane  corresponding  to  its  color.  This  permutation  of  points  in  a  column  ensures 
that  sources  and  destinations  end  up  in  the  same  plane  and  that  at  most  one  source 
and  one  destination  get  assigned  to  each  point  in  a  given  plane.  □ 

Lemma  3.2.1  now  follows  as  a  corollary  of  Lemma  3.2.2: 

Proof  of  Lemma  3.2.1:  The  idea  here  is  to  route  all  of  the  source-destination 
pairs  into  a  sub-cube  of  G  of  size  [mi]d,  and  then  apply  Lemma  3.2.2.  Hence  it 
suffices  to  show  that  the  sources  and  destination  can  be  packed  into  such  a  sub-cube 
using  each  edge  of  G  at  most  2  [mi]  times.  Again  pick  a  dimension  and  consider 
the  nd~l  planes  connected  by  columns  along  that  dimension.  Since  there  are  only 
m  point  to  go  around,  at  most  mi  planes  contain  more  than  [mi]J_1  points.  Call 
a  plane  with  more  than  that  number  heavy.  Clearly,  each  heavy  plane  contains  a 
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point  in  a  column  which  is  at  most  fma]  away  from  a  light  plane  with  an  empty 
space  in  that  column.  Route  this  point  along  its  column  to  that  empty  space  and 
proceed  in  this  manner  until  there  are  no  heavy  planes.  Since  all  routing  occurs  along 
columns,  and  each  column  intersects  at  most  [ma]  heavy  planes,  no  edge  gets  used 
more  than  fma' )  times.  Now  consider  the  sub-cube  of  size  [ma"]1*  chopped  up  into 
n  non-overlapping  slices,  corresponding  to  the  points  in  each  of  the  n  planes.  By 
induction  we  can  pack  the  points  of  each  plane  into  the  right  positions  indicated  by 
the  plane’s  corresponding  slice  using  each  edge  of  the  plane  at  most 
=  2[m^]  times.  Then  squash  down  along  the  columns  to  pack  all  of  the  slices  into 
the  sub-cube.  We  have  now  used  each  edge  in  the  column  dimension  at  most  2[m^] 
times.  □ 

By  letting  d  vary  from  1  to  log  N,  and  observing  that  the  edges  of  the  d- 
dimensional  grid  can  be  divided  into  2d  matchings,  Lemma  3.2.1  and  the  proof  of 
Theorem  1.2.1  indicate  that  any  embedding  of  a  d- dimensional  grid  in  the  plane  is 
within  an  O(diVa)  factor  of  optimal. 


3.3  Finding  Dense  Subsets  of  Points  in  the  Plane 

Let  fopt  :  G  — ♦  P  be  an  optimal  cost  embedding  of  a  d-dimensional  grid  among  the 
N  points  in  the  plane.  We  wish  to  prove  that  given  any  K  point  subset  S  C  P,  at 
least  half  the  points  of  5  are  contained  in  a  dense  sub- region  of  the  plane.  For  our 
purposes  dense  will  be  defined  in  terms  of  fopt(G)  and  K  by  the  following  lemma: 


Lemma  3.3.1  Given  any  fixed  K  point  subset  of  the  N  points  in  the  plane,  there 
exists  a  rectangular  sub-region  of  the  plane  containing  y  of  these  points  with  longest 
side  length  R,  where  R  <  0  _ 


Proof:  Consider  the  s  x  s  square  region  containing  all  of  the  points.  Starting  at  the 
left  side  of  this  region,  move  a  vertical  line  to  the  right  until  exactly  y  among  the  K 
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Figure  3-2:  Finding  dense  subsets  of  points  in  the  plane. 

distinguished  points  are  contained  in  the  left  region  of  the  square  defined  by  the  line. 
Do  the  same  with  a  vertical  line  working  from  the  right  to  the  left.  Now  let  R  be  the 
distance  between  these  two  vertical  lines.  See  Figure  3-2.  Label  the  points  in  the  left 
section  from  1  to  y .  Do  the  same  for  the  right  section.  By  Lemma  3.2.1  we  can  draw 
paths  in  the  optimal  d-dimensional  grid  connecting  corresponding  points  so  that  no 
edge  of  the  optimal  grid  is  used  more  than  O(K^)  times.  Each  of  these  paths  must 
cross  distance  R.  Hence  ||/opt(G)||  >  U  and  therefore  R<  O  Now 

do  the  same  procedure  with  horizontal  lines  from  the  top  and  bottom.  The  center 
rectangle  now  has  at  least  y  points  and  we  may  shrink  it  slightly  until  it  holds  exactly 

f  □ 

Our  embedding  procedure  will  be  as  follows.  Use  Lemma  3.3.1  to  isolate  y  points 
in  a  dense  region  5j.  Then  embed  an  y-node  piece  of  the  grid  G\  in  these  points.  Of 
the  remaining  points,  use  Lemma  3.3.1  to  isolate  y  in  another  perhaps  slightly  less 
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O 

o 

o 


o 

o 

o 


Figure  3-3:  G{  embeds  in  the  points  contained  in  5,  \  5,_ 

dense,  region  52.  Then  embed  an  ^-node  piece  of  the  grid  G2  in  these  points.  Proceed 
in  this  manner  until  all  of  the  points  in  the  plane  are  used  up.  Since  each  rectangular 
region  5,  is  obtained  by  sliding  vertical  and  horizontal  lines  in  from  the  boundary  of 
the  plane  until  the  right  number  of  points  are  captured,  Si  C  S2  C  •  •  •  C  Sjogw+i-  G, 
will  be  embedded  in  the  points  contained  in  S’,  \  5,_j.  See  Figure  3-3. 

In  this  manner,  we  need  a  strategy  for  dividing  G  up  into  log  TV  -f  1  pieces  of 
geometrically  decreasing  size.  Let  G  =  G\  U  G2  U  •  •  •  U  Giogyv+i  where  G\  is  obtained 
by  chopping  G  in  half  along  one  dimension,  and  each  G,  is  obtained  by  chopping  G,_i 
in  half  along  its  longest  dimension.  We  call  this  the  geometric  decomposition  of  G. 
Figure  3-4  shows  the  geometric  decomposition  of  the  two  dimensional  grid.  It  is  not 
hard  to  see  that  in  d  dimensions,  this  decomposition  has  the  following  property: 

Observation  3.3.2  Let  G  be  a  d-dimensional  grid  on  N  points.  Then  the  geometric 
decomposition  ofG  =  GiUG2U-  •  •UGjogyv+i  has  the  property  that  at  most  — 
edges  of  G  connect  G ,  to  Gj  for  j  >  i  and  for  j  >  i  +  d  no  edges  connect  G,  and  G;. 

Proof:  This  proof  requires  that  we  give  a  formal  definition  of  the  geometric  decom¬ 
position  of  the  grid  in  d  dimensions.  Label  the  pieces  in  the  geometric  decomposition 
G  =  Gi  U  G2  U  •  •  ■  U  GiogN+i-  Imagine  that  the  d-dimensional  grid  composed  of  the 
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Figure  3-4:  A  geometric  decomposition  of  the  grid, 
integer  coordinate  points  in  euclidean  d-space: 


G  —  {1,2, . . . ,  n}  x  •••  x  {1,2, ...  ,n} 

S .  v  •  - . . 

d 

Edges  connect  adjacent  points  in  the  grid  whose  coordinates  differ  along  only  one 
dimension.  Then  we  have: 


Gi  =  [1,^]  x  [l,n]  x  •••  x  [l,n] 

G2  =  [|  +  l,n]  x  [1,^]  x  [l,n]  x  ...  x  [l,n] 


r  ti  _  rTi  _  r  n . 

Gd  =  [g  +  !*«]  *  *  [g  +  hn]  x  [1,  -] 

Gd+i  =  [g  +  i,  — ]  x  [-  +  l,n]  x  •••  x  [-  +  l,n] 

/o  r3n  .  .n  3n,  ,n  ,  ,  ,n 

G*+  2  =  [—  +  Tn]  x  [-  +  i,  — ]  x  [-  +  l,n]  x  •-  x  [-  +  l,n] 


2^J+1  -  l  r2Lil+1  -  1 

° '  =  t-liijV.  "  +  1’")x---xl  n  +  1'n) 


imodd 
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x 

x 


2W  -  1 

2W  _  l 
1  2L3J 

V  -  -  -- 


n  +  1, 


2UJ+1  -  1 
2L3J+1 


■n  +  1,  ra]  x  •  •  •  x 


n] 

f2«J  -  1 
1  2liJ 


n  +  1,  n) 


(d-l)-(imodd) 

The  edges  of  G  which  leave  G,  and  go  into  smaller  Gj  for  j  >  i  all  connect  points 
along  the  ((i  —  1)  mod  d)  +  1  dimension.  Notice  that  for  each  progressively  smaller 
piece,  G,+i,  6r,+2,  •  •  •>  the  size  of  the  face  adjacent  to  G,  is  halved.  Since  (^r)^  is 
an  upper  bound  on  the  number  of  edges  leaving  G,  and  entering  G,+i,  the  first  part 
of  Observation  3.3.2  follows.  Furthermore,  if  we  look  at  the  ranges  of  coordinate 
values  in  the  formula  for  Gi,  we  see  that  by  level  G,+d,  all  of  these  ranges  have  been 
reduced  in  half,  and  hence  the  edges  leaving  Gi  cannot  possibly  enter  any  of  the  Gj 
for  j  >  i  +  d.  □ 

We  now  indicate  how  to  embed  Gi  into  Si  \  5,_i  using  a  divide  and  conquer 
subroutine.  For  simplicity’s  sake,  we  suppose  that  Gi  is  an  Af-node  d- dimensional 
grid.  It  will  be  clear  how  to  modify  the  subroutine  to  handle  the  case  where  G,  is  not 
perfectly  square. 


3.3.1  The  Divide  and  Conquer  Subroutine 

Let  Rj  be  the  length  of  the  longest  side  of  the  region  5,-  containing  the  Af-point  subset 
Pi  C  P  into  which  we  are  embedding  G,.  A  straightforward  method  of  embedding 
the  d-dimensional  grid  Gi  into  P,  consists  of  dividing  the  points  in  half,  recursively 
embedding  half  of  G,  on  each  set  of  points,  and  then  drawing  in  the  connecting 

edges.  If  we  axe  careful  about  how  we  divide  up  Pi,  we  can  achieve  an  approximation 
algorithm  which  yields  an  embedding  with  cost  bounded  by  a  function  of  M  and 
Ri.  In  order  to  implement  such  a  divide  and  conquer  algorithm,  we  use  the  natural 
decomposition  tree  for  the  grid.  In  two  dimensions,  we  simply  split  G,  down  the  middle 
to  get  two  rectangles.  Then  split  these  rectangle  down  the  middle  to  get  two  squares, 
proceeding  in  this  manner  for  logM  levels.  See  Figure  3-5.  This  idea  generalizes 
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I - 1  I - II - 1  □□  □□ 

| - 1| - 1  □□□□ 

I _ I  I _ II _ I  □□  □□ 

Figure  3-5:  The  natural  decomposition  tree  for  the  grid  is  used  in  a  divide  and  conquer 
embedding. 

to  d  dimensions  in  the  obvious  fashion.  The  important  observation  concerning  this 
decomposition  is: 

Observation  3.3.3  Two  adjacent  blocks  on  the  i*h  level  of  the  natural  decomposition 
tree  of  G  have  at  most  O  edges  of  G{  passing  between  them. 

In  order  to  prove  Observation  3.3.3,  we  give  a  formal  definition  of  the  natural  de¬ 
composition  tree  in  d  dimensions.  Let  G  be  an  N  node  d-dimensional  grid.  We  define 
the  standard  binary  decomposition  tree  on  G  by  recursively  splitting  subsections  of 
G  in  half.  For  example,  if  N  =  nd  and  G  is  the  d-dimensional  grid  composed  of  the 
integer  coordinate  points  in  euclidean  d-space: 


G  =  {l,2,...,n}  x  •••  x  {l,2,...,n} 

v  V  ' 

d 

At  the  first  level  we  split  G  into  two  rectangular  sub-grids  isomorphic  to 

{1,2,...,^}  x  {1,2,. . . ,n}  x  •  •  •  x  {1,2, ...  ,n} 

- - - - - - - 

d 

In  this  manner,  at  the  ith  level,  G  has  been  decomposed  into  2*  rectangular  grids,  all 
isomorphic  to 


{i,2, . . . ,  2^1 }  x  •••  x  t1’2’---’  2fil*  x  *  *■'  * 


t  mod  <4  d— t  mod  d 

Proof  of  Observation  3.3.3:  At  level  i,  where  blocks  have  size  we  split 
them  by  cutting  along  the  d  —  1*‘  dimensional  plane  perpendicular  to  the  t  mod  d-(- 1 
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We  are  almost  ready  to  prove  a  bound  on  the  cost  of  am  embedding  given  by  the 
divide  and  conquer  algorithm.  First,  however,  we  need  the  following  technical  lemma: 

Lemma  3.3.4  Given  N  points  in  an  s  x  s  square  region  of  the  plane,  for  0  <  i  < 
log  N,  we  can  subdivide  the  plane  into  disjoint  rectangular  regions  B\, ...,  B'2,,  each 
containing  £  points  so  that  perimeter  (B'-)  <  4\/2Vs.  Furthermore,  Z?‘  =  Bjj-iU 

Proof:  To  begin,  assume  the  plane  has  been  rotated  so  that  no  two  points  share 
the  same  vertical  or  horizontal  coordinates.  We  now  proceed  by  induction.  Let 
f?°  =  the  s  x  s  square  containing  all  the  points.  Then  pertmeter(B^)  =  4s. 
On  the  even  levels,  we  will  divide  boxes  by  drawing  horizontal  lines,  and  on  the 
odd  levels,  by  drawing  vertical  lines.  So  suppose  that  Z?2', ...,  B%1,  have  been  con¬ 
structed  to  satisfy  the  lemma.  Then  we  divide  Bj'  in  half  by  drawing  a  verti¬ 
cal  line  which  splits  the  points  contained  in  Bf  into  two  equal  sets  contained  in 
boxes  B\'£\  an<l  Bg+1.  Repeat  this  process  on  these  two  boxes,  using  horizon¬ 
tal  lines,  in  order  to  generate  B2^,  B24f+2l),  B2^,  and  B24f+1).  See  Figure  3- 
6.  We  then  have:  £^j+,)  perimeter  (Bf,+^)  =  £/LoPer*me*eK^4j’-^)  - 
T£jf2l2perimeter{B*')  <  2(4\/22‘s)  =  4\Z22<’+1)s.  The  lemma  follows  by  induction. 
□ 
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Figure  3-6:  Dividing  Bk  2  into  four  rectangles  containing  an  equal  number  of  points 

The  following  lemma  shows  that  this  divide  and  conquer  subroutine  gives  a  bound 
on  the  cost  of  embedding  Gx  into  P,  in  terms  of  M  and  Rt. 

Lemma  3.3.5  Suppose  we  are  given  M  points  in  region  Si  with  longest  side  length 
Ri.  Then  in  0(M log  M)  time  we  can  embed  the  d-dimensional  grid  G{  on  these 
points  costing  0(y/M log  MRT)  when  d  =  2  and  0(^2^ Ri)  when  d  >  2. 

Proof:  As  in  the  proof  of  Lemma  3.3.4,  divide  Si  into  regions  2?’,  0  <  i  <  logM, 
1  <  j  <2*.  Contained  within  the  rectangular  regions  Bk  at  level  j,  we  draw  in 
the  edges  from  the  (j  -(-  2)nd  level  of  the  natural  decomposition  tree  of  Gt  indicated 
in  Figure  3-5.  For  any  rectangular  region  Bk,  the  maximum  distance  between  any 
two  points  is  at  most  half  the  perimeter.  Hence  by  Observation  3.3.3,  drawing  in  all 
of  the  edges  at  level  j  will  cost  at  most  EL,  diam(Bi)  <  = 

|  2  VMRi  if  d=2 

“  2M^(2~^yRi  if  d>2 

Summing  over  the  log  M  levels  gives  us  the  lemma.  In  the  d  =  2  case,  the  sequence 
is  not  geometric  so  we  get  the  extra  log  factor.  We  leave  the  running  time  analysis 
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to  Section  3.4.  □ 

We  are  now  ready  to  give  the  main  algorithm  for  embedding  the  entire  d- 
dimensional  grid  G  into  P. 

3.3.2  The  Grid  Embedding  Algorithm 

1.  Use  Lemma  3.3.1  to  draw  log  TV  nested  rectangular  regions  Si  C  S2  C  •  •  •  C 
SiogN+i  in  the  plane  so  that  S\  contains  y  points  of  P,  Si  \  5,_ i  contains  ~ 
points  of  P,  and  diameter(Si)  =  0(||/Opt(G)||(|)-)^’')-  See  Figure  3-3. 

2.  For  i  =  1,  •  •  • ,  log  JV  +  1  use  the  divide  and  conquer  subroutine  to  embed  piece 
G{  from  the  geometric  decomposition  into  5,  \  5<_j. 

This  algorithm  gives  an  assignment  of  the  nodes  of  G  to  the  points  of  P.  Call  this 
map  h  :  G  —*  P.  The  bounds  on  the  diameters  of  the  nested  rectangles  S\  C  S2  C 
•  •  •  C  SiogN+i  allow  us  to  prove  the  following  lemma. 

Lemma  3.3.6  Given  N  points  in  the  plane  and  G  =  Gj  U  G2  U  •  •  •  U  Giogjv+i>  a 
geometric  decomposition  of  the  d-dimensional  grid,  we  can  embed  the  pieces  G\,  ..., 
Giog/v+i  in  the  plane  so  that  each  Gi  lies  inside  Si  and  the  total  cost  of  these  embed¬ 
dings  is  0(\\fopt(G)\\  log2  N)  ford  =  2  and  0(||/opt(G)||  log  N)  for  d>  2. 

Proof:  Draw  the  nested  rectangles  as  in  Figure  3-3.  Now  apply  the  divide  and  con¬ 
quer  algorithm  to  embed  Gi  on  the  points  in  S,  \  5,-_ j.  By  Lemma  3.3.5,  the  cost  of 
embedding  Gi  is  diameter(Si))  when  d  >  2  and  log  -^diameter (Si)) 

when  d  =  2.  Since  diameter(Si)  =  0( jj/opl(C7)|j(-^y)^i )  we  have  ||/i(Grt)||  = 
0(\\Ut(G)\\)  if  d  >  2  and  ||fc(G.-)||  =  (||/0Ft(G)||  log  N)  if  d  =  2.  Summing  over 
the  log  TV  levels  yields  Lemma  3.3.6.  □ 


All  that  now  remains  for  us  to  do  is  account  for  the  cost  of  the  edges  hooking  together 
the  h(G,)’s  and  we  will  have  an  upper  bound  on  the  cost  of  the  embedding  h(G).  We 
put  in  these  connecting  edges  in  log  N  steps.  The  ith  step  consists  of  putting  in  all 
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edges  which  extend  from  G,  to  Gj  where  j  >  i.  These  edges  are  contained  in  and  by 
Observation  3.3.2  there  are  0(^rr(^-)i31 )  of  them.  The  cost  of  connecting  G,  to  G:  is 
then  0  )  =  0  fe^H|/„,,(G)||)  =  0  (2'‘-at1»'-'»|!/„p,(G)||). 

For  j  =  i,-  •  • ,  i  +  d  this  sequence  is  geometrically  decreasing  when  d  >  2.  Hence,  the 
total  cost  of  the  ith  level  interconnecting  edges  is  0( ||/opt(G)j|).  Summing  over  all  i 
levels,  we  see  that  these  interconnect  edges  add  at  most  0(||/opt(G)||  log  N)  to  the 
total  cost.  This  finishes  the  proof  of  the  most  fundamental  result  in  Part  I:  Theorem 
1.2.2. 


3.4  Running  Time  Analysis 

We  can  implement  an  algorithm  to  find  the  nested  rectangles  by  drawing  vertical  and 
horizontal  lines  to  remove  the  required  fraction  of  points  from  the  margins  and  then 
incrementally  shrinking  the  size  of  the  resulting  rectangle  to  capture  the  exact  number 
of  points  needed.  This  will  require  us  to  sort  the  points  once  by  x  coordinate  and 
once  by  y  coordinate.  After  sorting,  each  region  S,  can  be  found  in  0(|r)  time.  Hence 
we  can  compute  S\  C  S2  C  ■  •  •  C  S\ogN+i  in  G(;VlogiV)  time.  Computing  the  map 
h  :  G  — *  P  is  accomplished  by  calling  the  divide  and  conquer  subroutine  to  map  each 
G,  int  ■>  5,.  Let  T(Mi)  be  the  running  time  for  this  subroutine  on  the  3/,-node  giid 
G, .  If  we  start  by  sorting  the  Mt  points  by  x  coordinate  and  by  y  coordinate,  at  the 
jth  level  we  can  then  split  each  region  B3k  in  0(jf)  time.  Hence,  if  we  assume  the  M, 
points  are  first  properly  sorted:  T(M,)  =  0{M,)  4-  4 T(^)  =  0(M,  log  A/,).  Adding 
in  the  initial  sorting  time,  we  still  have  T(Mt)  =  0(3/,  log  A/,).  Then  summing  over 
all  Gi ,  where  A/,  =  we  see  that  h  can  be  computed  in  0(A^  log.V)  time. 
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Chapter  4 


Arbitrary  Weighted  Graphs 


We  now  turn  to  the  case  where  we  wish  to  embed  an  arbitrary  weighted  graph  G  in 
Rd.  As  indicated,  we  are  able  to  give  polynomial  time  approximation  algorithms  for 
this  problem  in  the  cases  where  P  is  an  array  of  points  (Theorem  1.2.3)  or  a  randomly 
chosen  uniform  distribution  of  points  (Theorem  1.2.4). 

4.1  Embedding  into  an  Array  of  Points 

This  section  provides  the  proof  of  Theorem  1.2.3. 

Theorem  1.2.3  There  exists  a  polynomial  time  algorithm  which,  when  given  an 
arbitrary  N -node  weighted  graph  G  and  a  d-dimensional  array  of  N  points  P  in  Rd, 
embeds  G  in  P  with  cost  within  an  0(log2  N)  factor  of  optimal. 

In  this  section,  we  assume  that  the  points  P  are  arranged  in  an  evenly  spaced 
d-dimensional  square  array  in  Rd.  Let  G  =  (V,E)  be  an  arbitrary  weighted  N- 
node  graph  with  edge  weights  w(e)  for  each  e  €  E.  A  £  -  separator  of  G  is  a 
decomposition  of  V  into  disjoint  sets  U  and  W  such  that  min(\\U\\,\\W\\)  >  y.  Let 
S  be  the  set  of  edges  that  join  points  in  U  to  points  in  W.  Then  we  define  the  cost 
of  the  separator  to  be  ||5j|  =  5ZegsU>(e).  Leighton  and  Rao  [38]  have  recently  given 
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an  approximation  algorithm  for  finding  minimum  cost  separators.  We  will  use  the 
following  result  from  their  paper. 

Theorem  4.1.1  (Leighton,  Rao)  There  exists  a  polynomial  time  algorithm  SEP 
such  that  given  a  weighted,  undirected  graph  G,  SEP  finds  a  -F  —  -F  separator  for  G 
which  has  cost  less  than  0(\ogN)  times  the  optimal  cost  g  —  |  separator  of  G. 

We  are  now  in  a  position  to  prove  Theorem  1.2.3.  We  give  the  proof  for  the  case 
when  P  is  a  2-dimensional  grid  of  points.  It  is  not  hard  to  extend  the  ideas  to  handle 
the  general  case  where  the  dimension  d  is  a  constant  greater  than  2.  Let  OPT  be  the 
set  of  edges  crossing  the  optimal  1  |  separator. 

We  proceed  with  a  divide  and  conquer  assignment  of  points  in  G  to  the  square 
array  P.  Consider  the  optimal  embedding  fopt  :  G  —*  P.  Let  s  be  the  side  length  of 
the  array  P. 

Claim  4.1.2  s  =  0(^0^) 

Proof:  Move  a  line  perpendicular  to  one  dimension  of  P  from  left  to  right  across  the 
plane,  until  it  splits  off  the  leftmost  |  of  the  nodes  in  P.  From  this  starting  point, 
continue  to  move  the  line  across  the  array  until  it  splits  off  |  of  the  nodes.  At  each 
point  in  its  traversal,  the  line  cuts  >  j|OPT|j  edge  weight  from  fopt{G).  Since  the  line 
moves  distance  ^ ,  \\fopt(G)\\  >  ^\\OPT\\.  □ 

Now  use  SEP  from  Theorem  4.1.1  to  obtain  a  ^  —  75  separator  V  =  U  U  IF  of  G 
which  has  cost  within  O(log  N)  times  ||OPr||.  Let  C  be  the  set  of  edges  joining  U 
and  W.  We  can  now  write  G  =  (U  U  IF,  E\  U  E2  U  C).  Again,  move  a  line  from  left  to 
right  across  P  and  stop  at  the  point  where  the  line  cuts  off  exactly  ||f/||  >  £3  leftmost 
points.  (The  line  may  contain  a  jog.)  See  Figure  4-1.  Recursively  embed  (U,Ei)  and 
(IF,  E2)  in  the  left  and  right  sub- rectangles  of  P  defined  by  the  line.  Notice  that  since 
at  each  level  we  are  splitting  sub-graphs  of  G  with  F  _  2L  separators,  as  long  as  we 
cut  the  sub-rectangles  of  P  along  their  longest  side,  the  aspect  ratio  of  any  rectangle 
never  exceeds  10  to  1.  We  measure  the  cost  of  this  recursively  defined  embedding  by 
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Figure  4-1:  Embedding  G  into  a  grid  of  points. 

adding  up  the  cost  of  the  edges  crossing  cuts  at  each  level.  At  the  top  level,  the  cost 
of  embedding  C  in  P  is  bounded  above  by  ||C||s  <  |j/opt(G)||  log  iV.  To  bound  the 
costs  introduced  at  the  lower  levels,  we  will  need  a  slightly  more  general  claim. 


Claim  4.1.3  Let  s  be  the  length  of  the  longest  side  of  an  M-node  sub-rectangle  R 
of  P  with  aspect  ratio  bounded  by  some  constant  k.  Then  if  ( U,EX )  is  any  M-node 
subgraph  of  G  with  optimal  |  —  j  separator  O  PT ,  s  =  Q(^jfopfjj^)- 

Proof:  Move  a  line  from  left  to  right  across  P,  starting  where  the  line  first  splits  off 
Y  of  the  points  of  fopt{U)  and  stopping  when  it  splits  off  the  leftmost  ~  points.  Let 
the  distance  between  these  two  cuts  be  /j.  Repeat  this  process  with  a  horizontal  line 
moving  from  the  top  to  bottom  of  P,  and  let  the  distance  measured  this  way  be  /2- 
Then  >  ~  of  the  points  of  fopt{U)  are  contained  in  an  lx  x  l2  region.  Hence  since  R  is 
roughly  square,  s  =  0(max(li,  l2)).  Now  by  arguments  similar  to  those  given  above, 
\\Ut(U)\\  >  li\\OPT\\  for  i  =  1,2.  It  follows  that  s  =  Q(Ji jjgffi}11).  □ 

To  finish  the  proof  of  Theorem  1.2.3  consider  the  ith  level  of  recursion  in  our 
embedding  algorithm.  Let  P  =  P,  U  •  •  •  U  Pk  be  the  disjoint  union  of  the  points 
contained  in  the  k  =  2'  rectangular  regions  defined  at  this  level.  Let  P,  have  longest 
side  length  S{.  Then  let  V  =  U\  U  •  •  •  U  Uk  be  the  disjoint  union  such  that  U,  is 
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embedded  in  P,.  Let  OPT,  be  the  optimal  |  ^  separator  of  the  subgraph  of  G 

determined  by  l/;,  and  let  C,  be  the  edges  crossing  the  ^  separator  found  by 
SEP .  Then  the  cost  of  the  embedding  introduced  at  level  i  is  bounded  above  by 
<  Hi=i  Si\\0PTi\\0 (log  N).  By  Claim  4.1.3  this  sum  is  bounded  above 
by  £f-i  ||/op*(£A)||0(log  N)  <  \\fopt(G)\\0{\og  N).  Summing  over  all  log  N  levels  of 
the  recursion,  the  total  cost  of  the  embedding  is  jj/op2(G)||0(log2  N). 


4.2  Embedding  into  a  Randomly  Distributed  Set 
of  Points 

This  section  presents  the  proof  of  Theorem  1.2.4. 

Theorem  1.2.4  There  exists  a  polynomial  time  algorithm  which,  when  given  an 
arbitrary  N -node  graph  G  and  a  random  set  of  N  uniformly  distributed  points  P 
in  Rd ,  with  high  probability  embeds  G  in  P  with  cost  within  an  0(log2  N)  factor  of 
optimal. 

An  algorithm  similar  to  the  one  described  above  achieves  an  0(log2Ar)  times 
optimal  embedding  with  high  probability  when  the  points  P  are  randomly  distributed 
in  a  square  s  x  s  planar  region  5.  Again,  we  use  SEP  to  find  a  good  E  —  2b  separator 
of  G ,  split  the  planar  region  into  appropriately  sized  pieces  with  a  vertical  line,  and 
proceed  recursively  to  embed  the  two  sub-graphs  of  G  in  these  pieces.  It  is  not  hard 
to  show  that  with  high  probability  Claim  4.1.2  still  holds  in  the  randomly  distributed 
case.  Hence,  the  embedding  cost  generated  by  the  cut  at  the  top  level  of  the  algorithm 
with  high  probability  is  0(||/opt(G)||  log  N). 

All  that  is  needed  at  this  point  is  a  probabilistic  analog  of  Claim  4.1.3.  So  suppose 
we  are  at  the  ith  level  of  recursion,  and  we  are  embedding  some  M- node  subgraph 
(U,E\)  of  G  into  a  sub-rectangle  R  of  the  plane  with  longest  side  length  r.  As  in 
Claim  4.1.3,  consider  the  -y  points  of  fopt(Or )  contained  in  an  /j  x  l2  region.  It  still 
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holds  that  ||/opt(L’)||  >  /,|iOPT’||  for  i  =  1,2.  We  simply  need  to  show  that  with  high 
probability  r  =  0(max(li,  /2)). 

For  this,  we  will  need  to  prove  some  preliminary  probabilistic  results  about  uniform 
distributions  of  points  in  the  plane.  All  of  these  results  use  standard  Chernoff  bound 
analysis.  See  Section  7.5.1  for  a  review  of  Chernoff  bounds. 

We  call  a  sub-region  R  of  the  s  x  s  square  5  an  q  x  r2  sub-rectangle  if  the  sides 
of  R  are  parallel  to  the  sides  of  5  and  the  lengths  of  its  sides  are  rj  and  r2.  The 
expected  number  of  points  lying  in  R  is  Ep{R)  =  jV^r2 . 

Lemma  4.2.1  Fix  l  and  w.  Then  there  exists  a  constant  A  such  that  for  all  N  with 
high  probability  no  sub-rectangle  R  with  shortest  side  length  rj  such  that  w  <  r^  <  2 u\ 
longest  side  r2  such  l  <  r2  <  2/  and  area  >  x,2l°&N  contains  more  than  14EP{R)  or 
less  than  ■LEp(R)  points. 

The  proof  of  Lemma  4.2.1  will  use  the  following  claim: 

Claim  4.2.2  Fix  a  region  R  in  the  square  S  and  a  constant  c.  Then  there  exists 
a  constant  A  such  that  for  all  N  if  the  area  of  R  is  >  Xa2  -  then  with  probability 
>  1  —  N~c,  R  contains  <  | Ep(R)  or  >  lEp(R)  points. 

Proof:  Let  q  be  the  probability  that  a  point  lies  in  R.  Then  using  well  known 
Chernoff  bounds  [42]  we  know  that  Prob [P  contains  <  ^EP{R)  or  >  | EP(R )  points] 
<  e~0gN  for  some  constant  0.  Now  we  may  pick  A  sufficiently  large  so  that  q  >  -*°^  N 
where  a  is  large  enough  so  that  e-or/31o8/v  <  N~c.  □ 

Proof  of  Lemma  4.2.1:  The  lemma  is  vacuously  true  if  Iw  <  o{  tP*),  so 
we  may  assume  that  Iw  >  fl(*2l^jV)  and  <  0(j—^).  We  proceed  to  show  that 
very  few  rectangular  regions  R  are  either  so  dense  or  so  sparse  that  they  violate  the 

conditions  of  Lemma  4.2.1.  We  will  handle  the  case  where  R  contains  more  than 

2 

14Pp(P)  points  first.  Divide  S  into  sub-regions  of  size  l  x  w  in  two  ways.  One  way. 
the  length  l  side  is  horizontal,  and  the  other  way  it  is  vertical.  See  Figure  4-2  Then 
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Figure  4-2:  Divide  the  square  region  into  /  x  w  sub-rectangles  in  two  ways. 

2 

R  lies  inside  some  3x3  block  of  these  sub-regions.  There  are  less  than  2f^  of  these 
blocks.  So  the  probability  that  R  contains  >  14 EP(R)  points  is  <  2 
times  the  probability  that  some  fixed  3x3  block  contains  >  14 EP(R)  points.  Let 
B  be  some  such  fixed  block.  But  the  expected  number  of  points  for  such  a  block  is 
E'  =  9^N.  So  14 EP(R)  >  | E'.  But,  Prob[number  of  points  in  B  >  | E']  <  N~c  by 
Claim  4.2.2.  Hence  Prob[some  R  contains  >  14 EP(R)  points]  <  Ar-^c-1h 

Now  consider  the  case  when  R  contains  less  than  jgEp(R)  points.  Divide  S  into 
^  II  x  sub-regions.  Then  some  sub-region  lies  inside  of  R.  So  the  probability 
that  R  contains  <  j^Ep(R)  points  is  <  the  probability  that  some  one  of  these  small 
regions  S'  contains  <  ^  —rhT2,  points.  But  EP(S')  —  ,  hence  this  probability  is  < 

the  probability  that  S'  contains  less  than  ^ EP(S ').  Once  again,  by  the  above  claim, 
the  probability  that  a  fixed  S'  contains  <  ^EP(S')  is  <  N~c.  Hence  Probfsome  R 
contains  <  j^Ep{R)  points]  <  □ 

We  can  now  prove  the  stronger  Lemma  4.2.3  which  holds  for  any  sufficiently  large, 
bounded  aspect  ratio  sub- rectangle: 

Lemma  4.2.3  There  exist  a  constant  A  such  that  with  high  probability  no  sub- 
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rectangle  with  area  >  s  contains  more  than  14 EP(R)  or  less  than  Tj EP{R) 

points. 

Proof:  Let  R  be  a  sub-rectangle  satisfying  the  assumptions  of  the  lemma.  Let  r! 
and  r2  be  the  lengths  of  the  shortest  and  longest  sides  of  R.  Divide  the  ranges  of 
possible  values  for  rj  and  r2  up  into  sub-intervals  of  length  jj.  Each  pair  of  sub¬ 
intervals  corresponds  to  a  pair  of  values  for  l  and  w  in  Lemma  4.2.1.  There  are 
at  most  N 2  such  possible  pairs  of  values.  Notice  in  the  proof  above  that  A  does 
not  depend  on  /  or  w ,  so  we  may  pick  some  constant  A  to  suffice  for  all  and  yield 
probability  <  N~^c+2h  Since  there  are  <  N 2  sub-interval  pairs,  the  probability  that 
any  R  violates  the  conclusion  of  Lemma  4.2.3  is  <  N~c.  □ 

Lemma  4.2.4  With  high  probability,  any  sub-rectangle  R  generated  during  our  recur¬ 
sive  embedding  procedure  has  aspect  ratio  bounded  above  by  some  constant  k  whenever 
R  contains  >  A  log  N  points. 

Proof:  We  proceed  by  induction  on  the  construction  of  rectangles.  Lets  say  we  start 
with  some  rectangle  R'  which  has  aspect  ratio  <  k  and  contains  M  points  of  the 
distribution  within  its  area.  We  are  going  to  split  its  points  better  than  -E  —  with 
a  line  perpendicular  to  its  longest  side.  What  is  the  probability  that  -L  of  its  points 
are  squashed  up  into  some  sub-rectangle  R  with  aspect  ratio  >  k ?  See  Figure  4-3. 
When  such  a  squashing  occurs,  we  have  >  E  of  the  points  in  R'  packed  into  <  A 
of  the  area  of  R'.  Let  EV(R)  =  NAre^Rl  be  the  expected  number  of  points  in  the 
rectangle  R.  From  Lemma  4.2.3  we  know  that  with  high  probability  R'  has  area 
<  I877.S2.  This  implies  that  there  are  points  in  a  sub-rectangle  of  size  y^rs2.  For 
k  >  (18)(10)(14)  this  violates  Lemma  4.2.3.  □ 

We  now  return  to  our  original  goal:  proving  that  the  recursive  embedding  algo¬ 
rithm  using  SEP  generates  near  optimal  embeddings  with  high  probability  given  a 
uniform  distribution  of  points.  We  left  off  at  the  ith  level  of  recursion.  At  this  level, 
we  wish  to  embed  some  A/-node  subgraph  (U,  E\ )  of  G  into  a  sub-rectangle  R  of  the 
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Figure  4-3:  A  bad  distribution  of  points  leads  to  high  aspect  ratio  rectangles. 

plane  with  longest  side  length  r.  The  y  points  of  fopt(U)  are  contained  in  an  lx  x  l2 
region  R.  We  need  to  show  that  with  high  probability  r  =  0(mai(/1,/2)).  Now  if 
it  is  not  true  that  r  =  0(max(li,  I2)),  then  R  is  a  sparsely  populated  rectangular 
region  with  area  >  |r2  >>  lxl2  which  contains  only  M  points.  But  the  lxl2  region 
also  contains  >  y  points.  When  M  >  A  log  N  these  conditions  violate  Lemma  4.2.3. 

Now  we  must  handle  the  case  where  M  <  A  log  N.  Let  us  alter  the  algorithm  used 
for  embedding  arbitrary  graphs  in  arrays  so  that  when  regions  reach  size  <  A  log  N  we 
stop.  Then  we  have  sub-divided  V  =  Ux  U  •••  UUk  such  that  ||£/;||  <  AlogiV.  Further¬ 
more.  each  U,  is  contained  in  a  region  with  diameter  <  0  —  0(\\fopt{Ut)\\). 

Hence  we  can  embed  the  remaining  subgraphs  into  their  regions  any  way  we  want 
with  cost  <  ||/oPt(^)P2  log2  N  <  0(||/ope(G)||  log2  N).  Therefore,  with  high  prob¬ 
ability.  the  total  cost  of  the  embedding  is  0(||/opt(G)j|  log2  N).  □ 


Chapter  5 


Applications  to  Parallel 
Processing 


Several  important  practical  problems  in  parallel  processing  concern  the  efficient  em¬ 
bedding  of  large  computational  problems  into  parallel  architectures.  We  consider  a 
distributed  model  of  parallel  computing  where  a  parallel  algorithm  is  composed  of 
a  set  of  processes.  Each  process  performs  a  certain  set  of  computations  and  sends 
and  receives  data  from  other  processes.  To  avoid  confusion  of  processes  with  pro¬ 
cessors.  we  will  sometimes  refer  to  processes  as  tasks.  We  define  the  computation 
graph  for  an  algorithm  to  be  the  graph  which  has  one  node  for  each  process  and  edges 
between  nodes  whose  processes  exchange  data.  In  some  cases,  we  may  consider  the 
weighted  computation  graph  in  which  edge  weights  represent  the  relative  amount  of 
communication  carried  by  each  edge.  Examples  of  algorithms  and  their  corresponding 
computation  graphs  include: 

1.  Circuit  simulation.  Each  node  of  the  computation  represents  a  device  in  the 
circuit  with  edges  between  devices  which  exchange  electrical  information. 

2.  Numerical  methods  for  solving  differential  equations.  Various  finite  difference 
methods  yield  grids,  cubes,  and  other  computation  graphs  for  numerically  solv- 
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ing  differential  equations. 

3.  Fast  Fourier  transform.  The  natural  computation  graph  is  the  butterfly. 

A  problem  which  arises  naturally  in  this  context  concerns  how  to  assign  an  algorithm's 
tasks  to  a  machine’s  processors  and  efficiently  embed  the  computation  graphs  for  these 
algorithms  into  the  parallel  architecture.  Our  geometric  embedding  results  address 
this  problem  by  providing  algorithms  for  assigning  tasks  to  processors  for  a  large  class 
of  machines:  array- based  parallel  processors. 


5.1  Minimizing  Communication  Load 

We  measure  the  efficiency  of  an  embedding  in  a  parallel  architecture  in  terms  of 
the  communication  load  which  that  embedding  induces  on  the  interconnect  network. 
Communication  load  is  the  total  volume  of  inter-processor  communication  carried  by 
the  interconnect  network  on  the  architecture.  For  any  given  parallel  algorithm,  this 
load  may  vary  dramatically  depending  upon  how  tasks  are  assigned  to  processors. 
Let  mi, . . .  ,mk  be  the  messages  which  must  be  exchanged  between  processors  dur¬ 
ing  the  algorithm.  Then  if  d(m )  is  the  distance  traveled  by  a  message  through  the 
interconnection  network,  the  total  communication  load  is  simply:  In  this 

manner,  communication  load  is  a  measure  of  the  total  volume  of  traffic  the  network 
will  be  required  to  bear  during  the  algorithm. 

As  an  illustration,  consider  a  parallel  processing  problem  represented  as  a  graph 
on  a  set  of  tasks:  ffc.  The  weight  of  an  edge  in  this  graph  represents  the 

amount  of  communication  required  between  tasks  f,  and  tj.  See  Figure  5-1.  We 
need  to  optimally  assign  these  tasks  to  processors  in  a  manner  which  minimizes  the 
communication  load  on  the  network.  For  example,  if  we  are  interested  in  running  a 
circuit  simulation,  then  each  task  t{  would  correspond  to  a  device  in  the  circuit  and 
we  would  need  to  minimize  communication  load  induced  by  simulating  the  electrical 
signals  transmitted  between  devices. 
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Figure  5-1:  Computation  graph  for  a  problem  to  be  solved  on  a  parallel  machine. 

Our  results  cn  weighted  graph  embeddings  allow  us  to  give  a  provablv  near-optimal 
approximation  algorithm  for  minimizing  communication  load  on  array-based  parallel 
processors.  Here  we  understand  an  array-based  architecture  to  mean  any  iV-node  grid, 
cube,  or  higher  dimensional  array.  As  a  corollary  of  Theorem  1.2.3  we  can  embed  any 
weighted  computation  graph  in  an  array-based  machine  minimizing  communication 
load  to  within  an  0(log2  N)  factor  of  optimal. 


Corollary  5.1.1  There  exists  a  -polynomial  time  algorithm  which,  when  given  an  N- 
processor  array-based  machine  M  and  a  parallel  processing  task  consisting  of  a  set 
of  tasks  t\,...,tfi/  with  communication  costs  c(t,,t:)  between  tasks ,  can  compute  an 
assignment  of  tasks  to  the  processors  of  M  which  achieves  a  communication  load 
which  is  within  0(log2  N)  times  optimal. 


Proof:  In  a  fc-dimensional  array-based  processor,  the  distance  between  two  processors 
Pi  and  p2  is  dist(px,p2),  where  dist  is  the  Manhattan  metric  in  &-space.  Within 
constant  factors  depending  on  k.  this  metric  is  approximated  by  the  euclidean  distance 
function.  Let  G  =  (V,  E)  be  the  complete  graph  on  V  =  tx, . . . ,  tjv  with  edge  weights 
c(ti,  tj).  Then  we  simply  need  to  solve  the  weighted  geometric  embedding  problem  of 
mapping  G  into  a  ^-dimensional  grid.  □ 
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5.1.1  Simulating  Other  Architectures 

Corollary  5.1.1  can  also  be  interpreted  as  providing  an  algorithm  for  near-optimal 
simulations  of  other  architectures  on  array-based  machines.  An  alternate  architecture 
can  be  represented  as  an  iV-node  computation  graph  where  nodes  are  processors  and 
edge  weights  represent  the  amount  of  traffic  typically  present  on  the  link  between  two 
processors.  VVe  can  then  use  Corollary  5.1.1  to  embed  the  alternate  architecture  in 
the  array-based  machine  in  a  manner  which  minimizes  communication  load  to  within 
an  0(log2  N)  factor  of  optimal. 


5.2  Dynamic  Allocation  of  Resources  on  a  Mul¬ 
tiprocessor 

Suppose  we  are  given  an  Ar-processor  parallel  machine  configured  as  a  graph  G  — 
{V,  E)  where  J|V||  =  N,  and  a  schedule  of  problems  to  process:  pi, . . . ,  p*.  Each 
problem  p,  requires  a  certain  number  rii  of  processors.  The  machine’s  operating 
system  dynamically  allocates  processors  to  problems.  As  problems  finish,  holes  of 
idle  processors  open  up  in  the  network.  How  can  we  efficiently  assign  processors  from 
the  holes  to  the  remaining  problems  in  a  manner  which  minimizes  the  communication 
load  on  the  network?  If  each  problem  p,  corresponds  to  a  computation  graph  c,,  this 
problem  consists  of  finding  optimal  embeddings  of  c*  into  the  network  holes.  Using  our 
results  from  Theorem  1.2.4  on  embedding  weighted  graphs  into  randomly  distributed 
points  in  euclidean  space,  we  can  show  how  to  embed  an  arbitrary  computation 
graph  Conn  nodes  into  randomly  distributed  subsets  of  an  array  processor  within  a 
0(log2n)  factor  of  optimal. 

Corollary  5.2.1  There  exists  a  polynomial  time  algorithm  which,  when  given  a  k- 
dimensional  array  processor  M  with  n  randomly  distributed  idle  processors  and  a 
computation  graph  C  on  n  nodes,  with  high  probability  can  assign  nodes  of  C  to  the 
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idle  processors  in  a  manner  which  minimizes  the  additional  communication  load  on 
M  to  within  an  0(log2  n)  factor  of  optimal. 

Given  the  large  number  of  parallel  algorithms  devised  for  grids,  cubes,  hypercubes, 
butterflies,  and  shuffle  exchange  graphs,  we  are  frequently  interested  in  problems  p , 
which  have  one  of  these  graphs  as  their  computation  graph.  For  example,  algorithms 
for  solving  differential  equations  frequently  have  computation  graphs  which  are  grids, 
cubes,  or  higher  dimensional  arrays.  Fast  Fourier  transform  and  Batcher's  merge  sort 
run  naturally  on  the  butterfly.  In  these  cases,  we  can  drop  the  assumption  that  idle 
processors  are  randomly  distributed,  and  give  faster  algorithms  which  always  assign 
processors  in  a  near-optimal  manner.  In  fact,  when  C  is  a  hypercube,  shuffle-exchange 
graph,  or  butterfly,  any  processor  assignment  strategy  is  near-optimal. 

Corollary  5.2.2  Given  a  k-dimensional  array  processor  M  with  n  idle  processors 
and  a  computation  graph  C  on  n  nodes  which  is  a  uniformly  weighted  hypercube, 
butterfly,  or  shuffle  exchange  graph,  then  any  assignment  of  the  nodes  of  C  to  the 
idle  processors  creates  additional  communication  load  which  is  within  an  O(logn) 
factor  of  optimal. 

For  grids,  cubes,  and  higher  dimensional  arrays,  our  geometric  embedding  results 
yield  very  fast  approximation  algorithms  for  the  processor  assignment  problem. 

Corollary  5.2.3  There  exists  an  O(nlogn)  time  algorithm  vhich,  when  given  a  2- 
dimensional  array  processor  M  with  n  idle  processors  and  a  computation  graph  C 
on  n  nodes  which  is  an  uniformly  weighted,  can  assign  the  nodes  of  C  to  the  idle 
processors  to  create  additional  communication  load  which  is  within  an 

1 .  0(log2  n)  factor  of  optimal  when  C  is  a  grid. 

2.  O(logn)  factor  of  optimal  when  C  is  a  cube,  or  higher  dimensional  array. 
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5.2.1  Wafer-Scale  Integration  and  Reconfiguring  Around 
Faults  in  an  Array- Based  Processor 

Now  suppose  we  are  given  an  array-based  parallel  machine  with  a  number  of  faulty 
nodes.  Similarly,  consider  a  silicon  wafer  with  processors  arranged  in  an  intercon¬ 
nected  grid  where  some  of  the  processors  are  faulty.  We  assume  a  model  where 
messages  can  still  pass  through  a  failed  processor,  but  its  processing  capabilities  are 
dead.  Then  the  live  nodes  correspond  to  a  set  of  points  P  in  the  plane  or  Rd .  For  arbi¬ 
trary  faults,  Theorems  1.2.2  and  1.2.1  yield  algorithms  for  reconstructing  a  grid,  cube, 
hvpercube.  shuffle  exchange  graph,  or  butterfly  on  the  live  nodes  with  near-optimal 
communication  load.  On  the  other  hand,  if  the  faults  are  randomly  distributed,  then 
Theorem  1.2.4  yields  an  algorithm  for  embedding  arbitrary  architectures  on  the  live 
nodes  with  near-optimal  communication  load. 

5.3  Configuring  Large  Area  Distributed  Comput¬ 
ing  Networks 

Suppose  we  need  to  link  together  a  number  of  widely  dispersed  computing  centers 
into  a  single  network.  For  example,  the  weather  service  may  need  to  configure  a  large 
number  of  weather  stations  around  the  country  into  a  single  network  for  weather 
prediction.  We  may  wish  to  dynamically  add  and  remove  centers  from  this  network 
as  time  progresses.  For  example,  if  the  weather  service  were  tracking  and  analyzing  a 
storm,  the  stations  in  the  network  would  be  changing  over  time  as  the  storm  moved 
across  the  country.  In  the  case  of  such  large  area  distributed  networks,  communica¬ 
tion  costs  between  centers  will  comprise  a  significant  percentage  of  the  total  cost  of 
operating  the  network.  If  we  assume  that  communications  between  centers  will  take 
place  over  existing  telecommunications  channels,  the  operator  of  the  network  will  pay 
rates  proportional  to  the  distance  traveled  for  each  piece  of  communication.  Hence 
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the  operator  would  be  interested  in  configuring  the  network  so  as  to  minimize  the 
total  embedded  edge  length.  Assuming  that  centers  lie  on  the  plane,  this  reduces  to  a 
geometric  embedding  problem.  Our  results  indicate  how  to  produce  near-optimal  so¬ 
lutions  to  this  problem  when  the  network  is  a  grid,  cube,  hvpercube,  shuffle-exchange 
graph,  or  butterfly.  For  other  network  topologies,  we  can  achieve  the  same  results 
provided  the  centers  are  uniformly  distributed. 
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Part  II 

Query- Retrieval 
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Chapter  6 


Introduction 


6.1  Overview  of  Query- Retrieval  Results 

Query- retrieval  problems  have  received  considerable  attention  in  the  computational 
geometry  literature  [10,  13,  14,  15,  16,  18,  20,  24,  25].  Typically,  a  query- retrieval 
problem  consists  of  a  set  of  n  geometric  objects  (e.g.  points  in  the  plane),  and  an 
unlimited  equence  of  queries  (e.g.  “of  the  n  points  in  the  plane,  find  the  k  points 
that  are  closest  to  position  (x,  y)").  The  task  is  to  preprocess  the  n  objects  and  devise 
a  data  structure  so  that  each  query  can  be  answered  as  quickly  as  possible. 

The  most  important  measures  of  efficiency  in  a  query- retrieval  problem  are  the 
space  required  by  the  data  structure  and  the  time  needed  to  answer  each  query.  For 
example,  typical  query-retrieval  problems  require  fi(n)  space  and  Cl(k  +  logo)  time 
per  query  where  k  is  the  size  of  the  response  to  the  query.  The  algorithms  and  data 
structures  described  in  this  part  of  the  thesis  all  achieve  the  optimal  time  bound  and 
use  0(n)  or  0(nlogn)  space,  depending  on  the  problem  being  solved.  Preprocessing 
time  and  space  (used  to  construct  the  data  structure)  are  also  of  concern,  but  are 
usually  not  as  important  as  data  structure  space  and  query  response  time.  All  of  the 
algorithms  described  in  Part  II  of  this  thesis  use  polynomial  preprocessing  time  and 
space.  Using  probabilistic  methods,  the  preprocessing  time  and  space  can  be  reduced 
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to  0(n  log2  n). 

6.2  Main  Results 

In  this  part  of  the  thesis,  we  describe  a  new  technique  for  solving  query-retrieval 
problems  in  optimal  time  with  optimal  or  near-optimal  space.  The  technique  incor¬ 
porates  planar  separators,  filtering  search,  and  the  probabilistic  method  to  compact 
kth-order  Voronoi  diagrams  (and/or  other  suitable  proximity  diagrams)  from  k°(1)n 
space  to  O(n)  space  without  losing  any  of  the  information  that  is  essential  for  solving 
query  problems.  As  examples,  we  use  the  technique  to  construct  algorithms  and  data 
structures  for  A:-nearest  neighbor  search,  circular  range  search,  and  half-space  range 
search.  A  brief  description  of  each  of  these  problems  and  of  our  results  is  provided 
below. 


6.2.1  Planar  fc-Nearest  Neighbor  Search 

Input :  A  set  P  of  n  points  in  the  Euclidean  plane  E 2  and  an  integer  k  >  0. 

Query:  Find  the  k  points  in  P  that  are  closest  to  a  query  point  q. 

We  show  how  to  solve  this  problem  using  0(n)  space  and  0(k  -f  logn)  time  per 
query,  both  of  which  are  optimal.  If  k  is  not  fixed,  but  is  given  as  a  part  of  the 
query  (i.e.,  the  query  is  a  pair  (q,  k)  with  q  €  E2  and  k  <  n),  then  our  solution  uses 
O(nlogn)  space  and  O(fc  +  logn)  time  per  query.  The  best  previously  known  solution 
to  this  problem  is  due  to  Chazelle,  Cole,  Preparata,  and  Yap  [15]  who  construct  a 
data  structure  using  0(n(log  n  log  log  n)2)  space. 

We  also  consider  the  ^-nearest  neighbors  problem  in  some  other  non-euclidean 
metrics.  For  example,  we  show  how  to  construct  data  structures  and  algorithms  with 
equivalent  performance  in  the  power-distance  metric  and  the  additive-weight  metric. 
This  metrics  are  defined  and  described  in  Section  7.1.2.  Here  we  simply  point  out 
that  Voronoi  diagrams  in  these  alternative  metrics  are  well  known  and  have  several 
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applications  [6,  5.  7,  8,  22,  20.  41,  43],  but  the  A:-nearest  neighbor  searching  problem 
in  these  metrics  does  not  appear  to  have  been  studied  previously.  As  an  example 
of  the  applications  of  these  metrics,  we  use  our  technique  for  computing  ^-nearest 
neighbors  in  the  power-distance  metric  to  solve  the  3-dimensional  half-space  range 
search  problem. 

In  the  special  case  that  each  query  point  q  belongs  to  a  restricted  set  Q  of  0(n) 
points  (e.g.,  Q  =  P),  the  query  response  time  of  our  algorithms  can  be  improved  to 
O(k)  without  increasing  the  space.  This  restricted  searching  result  could  be  useful 
in  conjunction  with  heuristics  for  large  traveling  salesman  problems  that  repeatedly 
compute  the  k  nearest  neighbors  of  various  cities  along  the  tour  [34].  For  small 
problem  instances,  it  is  possible  to  simply  precompute  the  ^-nearest  neighbors  of 
every  node  and  store  them,  thereby  using  0(kn )  space.  For  large  n,  however,  this 
approach  is  not  possible,  and  the  techniques  described  in  Part  II  of  this  thesis  become 
more  appropriate.  In  this  application,  our  approach  uses  O(n)  space  and  O(k)  query 
time,  both  of  which  are  optimal. 

6.2.2  Circular  Range  Search 

Input:  A  set  P  of  n  points  in  the  Euclidean  plane  E2. 

Query:  Find  all  points  of  P  contained  in  a  disk  in  E 2  with  radius 
r  centered  at  q. 

We  show  how  to  solve  this  problem  using  0(n  log  n)  space  and  0(k+log  n )  time  per 
query  where  k  is  the  number  of  points  in  the  disk.  Several  algorithms  for  this  problem 
have  been  reported  in  the  literature  [10,  13,  15,  16].  The  best  previously  known 
algorithm  (due  to  Chazelle,  Cole,  Preparata,  and  Yap  [15])  uses  0(n(logn  log  log  n)2) 
space  and  0(k  -f  lo gn)  time  per  query.  We  transform  the  circular  range  search 
problem  into  the  planar  ^-nearest  neighbors  search  problem  and  apply  the  filtering 
search  technique  given  in  [13]  to  obtain  the  result. 
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If  the  centers  of  the  query  disks  are  known  ahead  of  time  and  if  the  number  of 
centers  is  0(n)  then  the  query  time  can  be  improved  to  0(k)  without  increasing  the 
space. 

6.2.3  Half-Space  Range  Search  in  Three  Dimensions 

Input:  A  set  P  of  n  points  in  Euclidean  space  E 3. 

Query.  Find  all  points  of  P  that  are  contained  in  the  half  space  defined 
by  ax  +  by  +  cz  <  d. 

We  show  how  to  solve  this  problem  with  0(n  log  n)  space  and  0(k  +  log  n)  query 
time,  where  k  is  the  number  of  points  found  to  be  contained  in  the  half-space.  Previ¬ 
ously,  Chazelle  and  Preparata  [20]  had  given  a  data  structure  for  this  problem  using 
0(n  log5  n  log  log4  n)  space  which  was  later  improved  to  0(n  log2  n  log  log  n )  space  by 
Clarkson  and  Shor  [24].  We  transform  the  range  search  problem  into  the  k- nearest 
neighbor  search  problem  in  the  power-distance  metric  and  use  filtering  search  to  ob¬ 
tain  our  solution. 


6.3  Outline  of  Part  II 

The  remainder  of  Part  II  is  divided  into  three  sections.  Chapter  7  presents  our 
compaction  technique  by  constructing  an  0(n)-space  data  structure  for  fc-nearest 
neighbor  search  in  the  Euclidean  plane.  Initially  a  randomized  construction  is  given, 
and  later  (in  Section  7.S)  we  show  how  to  use  the  probabilistic  method  [42]  to  give 
a  deterministic  construction.  We  also  indicate  how  the  compaction  technique  can  be 
applied  to  nearest  neighbor  search  problems  in  other  metrics.  Section  7.10  briefly 
illustrates  how  we  can  build  a  linear  data  structure  for  ^-nearest  neighbor  search 
with  significantly  improved  preprocessing  time  by  re-introducing  some  randomness 
into  the  deterministic  construction.  Finally,  Chapter  8  discusses  applications  of  the 
compaction  technique  to  circular  range  search  problems  in  the  plane  and  half-space 
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range  search  problems  in  three  dimensions.  We  conclude  this  chapter  with  a  broad 
generalization  of  the  technique  and  suggest  applications  to  other  retrieval  problems. 
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Chapter  7 


Voronoi  Diagram  Compaction 

7.1  Review  of  Voronoi  Diagrams 

The  importance  of  Voronoi  Diagrams  in  computer  science  and  other  fields  has  long 
been  recognized  [1,  6,  36].  Here  we  give  a  brief  overview  of  the  properties  of  Voronoi 
diagrams  in  the  euclidean  plane  and  other  metrics  which  are  important  to  the  results 
proved  in  this  thesis.  Throughout  this  section  we  let  P  =  pi, . . .  ,pn  be  a  set  of  points 
in  the  plane  which  are  in  general  position. 

7.1.1  The  Euclidean  Plane 

For  any  pair  pi,pj  €  P,  let  h(pi,pj )  be  the  half- plane  containing  p,  which  is  defined 
by  the  perpendicular  bisector  between  p,  and  pj.  Then  the  Voronoi  diagram  defined 
by  pi, . . .  ,pn  is  the  partition  of  the  plane  into  convex  regions  R(pi), . . . ,  R(pn)  where 
R(p,)  =  D&i  h{Pji  Pi)-  1°  this  manner,  the  boundary  edge  between  any  two  convex 
regions  in  the  Voronoi  diagram  is  a  segment  of  a  perpendicular  bisector  between  two 
points  p,  and  p3.  R(pi)  consists  of  those  points  in  R?  which  are  closer  to  p,  than 
any  other  point  p:  €  P.  Likewise,  the  kth- order  Voronoi  diagram  for  pi,...,p„  is 
the  partition  of  the  plane  into  convex  regions  R(T)  where  T  =  p, pu  is  any 
subset  of  k  points  and  R(T)  =  f)PeT.qeP\T  h(p,q)-  R(T)  consists  of  all  points  in 
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the  plane  whose  k- nearest  neighbors  are  exactly  T .  Similarly,  the  boundary  edge 
between  any  two  convex  regions  in  the  kth -order  Voronoi  diagram  is  also  a  segment 
of  a  perpendicular  bisector  between  two  points  pi  and  pj.  This  yields  the  following 
observation: 

Observation  7.1.1  Whenever  R(T)  and  R(S)  are  adjacent  faces  in  the  kth-order 
Voronoi  diagram.  T  and  S  differ  in  exactly  one  point. 

Proof:  Let  s  and  t  be  the  points  whose  perpendicular  bisector  defines  the  common 
edge  between  R{T)  and  R(S).  WLOG  suppose  that  R{T)  C  h(t,s)  and  R{S)  C 
h(s,t).  Then  s  £  S  and  t  €  T  so  S  and  T  differ  in  at  least  one  point.  Now  suppose 
that  S  and  T  differ  on  another  pair  of  points  s'  and  t'.  Then  R{S)  C  /i(s,  t)  (T  h(s',  t') 
and  R(T)  C  h(t,  s)  f)  h(t',  s').  But  fi(s,t)f \h(s',t')  and  h(t,  s)  f)  h(t\  s')  do  not  share 
a  common  edge.  We  have  a  contradiction  and  S  and  T  must  differ  in  at  most  one 
point.  □ 

Along  with  the  above  observation,  this  thesis  uses  the  following  well  known  prop¬ 
erties  of  the  kth- order  Voronoi  diagram  V: 

1.  V  is  a  planar  graph  with  0{kn)  edges  and  faces.  [36] 

2.  V  can  be  computed  in  0(k2n  -f  nlogn)  time.  [1] 

7.1.2  Other  Metrics 

In  this  thesis,  we  also  apply  the  compaction  technique  to  solve  query-retrieval  prob¬ 
lems  in  other  metric  spaces.  In  particular  we  address  the  power-distance  metric  and 
the  additive-weight  metric.  Power-distance  and  additive-weight  metric  Voronoi  dia¬ 
grams  are  well  known  and  have  several  applications  [5,  6,  7,  8,  20,  22,  41,  43].  Below 
we  give  a  brief  introduction  to  the  properties  of  Voronoi  diagrams  in  these  metric 
spaces  which  are  important  to  the  results  appearing  in  this  thesis. 
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The  Power-Distance  Metric 

In  the  power-distance  metric ,  a  weight  wp  >  0  is  provided  with  each  point  p  G  P. 
The  distance  between  p  and  q  is  defined  to  be  pow(p.q)  =  d(p,q)2  —  w 2  where  d(p,q) 
is  the  euclidean  distance  between  p  and  q.  The  following  geometric  interpretation 
can  be  given  to  pow(p,q):  Let  L  be  the  line  through  q  that  is  tangent  to  a  disk  with 
radius  w2  and  center  p.  Then  \]pow(p,  q)  is  the  length  of  the  segment  between  q  and 
the  point  on  L  tangent  to  the  disk.  Given  two  points  p,q  6  P,  the  locus  of  points 
with  equal  weighted  distance  from  both  p  and  q  is  a  straight  line.  Hence,  the  concept 
of  kth- order  Voronoi  diagram  generalizes  nicely  to  this  metric.  Aurenhammer  [5,  6] 
and  Edelsbrunner  [27]  discuss  some  properties  of  kth- order  Voronoi  diagrams  in  the 
power-distance  metric,  and  they  call  these  diagrams  ktk- order  power  diagrams.  The 
properties  that  are  most  important  to  us  include  the  facts  that  the  edges  defining  the 
kth- order  power  diagram  are  straight  lines  and  that  the  diagram  contains  only  0(k2n) 
faces.  We  discuss  these  diagrams  in  more  detail  in  Section  8.2. 


Weighted  Metrics 

In  the  additive-weight  metric ,  the  distance  between  p  and  q  is  defined  to  be  d(p,q)+wp. 
The  ktfl- order  Voronoi  diagram  for  this  metric  (called  a  weighted  Voronoi  diagram)  has 
only  0(kn )  edges.  Rosenberger  [41]  and  Ash  and  Bolker  [4]  have  shown  that  such  dia¬ 
grams  can  be  computed  in  0(k2n  log  n)  time  using  0(kn )  space.  The  edges  bounding 
the  regions  of  kth- order  weighted  Voronoi  diagrams  are  second  degree  curves.  [41]  As 
it  turns  out,  we  are  able  to  apply  our  compaction  techniques  to  the  ^-nearest  neigh¬ 
bors  problem  in  the  additive-weight  metric.  This  serves  as  an  example  to  illustrate 
the  fact  that  the  techniques  presented  here,  and  fully  developed  for  the  euclidean 
metric,  can  be  extended  to  other  metric  spaces. 
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7.2  Planar  Point  Location 

In  the  construction  of  our  data  structure  for  query- retrieval  problems,  we  use  a  planar 
point  location  data  structure  due  to  Kirkpatrick.  [35]  The  planar  point  location  prob¬ 
lem  assumes  a  coordinate  embedding  of  a  planar  graph  (such  as  a  Voronoi  diagram! 
in  R2.  If  the  graph  has  n  edges  then  we  can  build  a  data  structure  which  can  answer 
the  following  query  in  O(logn)  time:  which  face  of  the  planar  graph  does  the  query 
point  lie  in  ?  Kirkpatrick’s  planar  point  location  data  structure  requires  only  O(n) 
space  and  can  be  constructed  with  <9(nlogn )  preprocessing  time. 

A  different  planar  point  location  technique  is  given  by  Edelsbrunner,  Guibas.  and 
Stolfi  [28].  This  technique  also  achieves  linear  space  and  logarithmic  time  planar 
point  location.  It  has  the  advantage  over  Kirkpatrick’s  technique  in  that  it  can  be 
generalized  to  to  solve  point  location  problems  in  sub-divisions  of  the  plane  that  are 
defined  by  curved  edges  such  as  hyperbolas.  The  fact  that  we  are  able  to  do  optimal 
time  and  space  planar  point  locations  in  sub-divisions  defined  by  curves  is  important 
in  latter  sections  of  this  thesis  where  we  extend  our  compaction  techniques  to  Voronoi 
diagrams  in  non-euclidean  metric  spaces. 

7.3  The  Voronoi  Diagram  Compaction  Technique 

At  this  point,  we  are  ready  to  present  the  main  results  of  Part  II  of  this  thesis. 
Let  V  be  the  triangulated  kth-OTdex  Voronoi  diagram  for  a  set  P  of  n  points  in  the 
plane  under  the  standard  Z,2  (euclidean)  metric.  We  assume  that  the  points  of  P  are 
in  general  position.  Then  V  is  a  planar  graph  with  0(kn )  edges  and  O(kn)  faces 
[36].  Using  Kirkpatrick’s  results  [35]  on  planar  point  location,  it  is  well  known  how  to 
preprocess  V  to  answer  ^-nearest  neighbor  queries  in  0(log  n  +  k)  time.  Unfortunately, 
the  resulting  data  structure  occupies  0(kn)  space.  In  this  section  we  construct  an 
optimal  data  structure  using  only  0(n )  space  which  achieves  the  optimal  0(logn  +  k) 
time  bound  for  answering  k- nearest  neighbor  queries. 
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7.4  Planar  Separator  Techniques 

We  begin  bv  taking  the  dual  of  V:  V.  Since  the  points  are  in  general  position.  V'  is 
a  planar,  degree  3  graph  with  O(fcn)  edges  and  nodes.  A  node  v  of  V  corresponds 
to  a  region  of  the  plane.  Any  query  point  in  this  region  has  the  same  set  Sv  of  k- 
nearest  neighbors  in  P.  It  follows  from  Observation  7.1.1  that  whenever  v  and  w  are 
adjacent  nodes  in  V,  Sv  and  Sw  differ  in  at  most  one  point.  For  each  edge  e  =  (t\  u •) 
in  V ,  label  e  with  the  set  of  points  Le  =  SVC\SW.  So  Le  contains  either  k  —  1  or  k 
points.  Now  fix  a  point  p  €  P  and  consider  the  faces  of  V  which  contain  p  in  their 
list  of  Avnearest  neighbors.  These  faces  correspond  to  a  set  of  nodes  in  V.  Let  (V)p 
consist  of  these  nodes,  together  with  the  edges  e  such  that  p  6  Le.  Then  we  have  the 
following  observation: 

Observation  7.4.1  For  all  p  6  P,  (V)p  is  connected. 

Proof:  Consider  a  node  v  €  (Vr)p  and  the  corresponding  face  Fv  in  V.  Draw  a 
straight  line  from  any  point  in  Fv  to  p.  Moving  along  this  line,  we  are  always  getting 
closer  to  p,  so  each  face  of  V  which  we  pass  through  contains  p  as  one  of  its  k- nearest 
neighbors.  Furthermore,  since  they  are  adjacent,  these  faces  which  the  line  passes 
through  correspond  to  a  connected  subset  of  nodes  in  ( V")p.  Let  w  be  the  node  in 
(C)p  whose  corresponding  face  in  V  contains  p.  Then  we  have  shown  that  v  is  in  the 
same  connected  component  as  w  for  all  v  £  {V)p.  □ 

At  this  point,  we  would  like  to  partition  V  into  groups  containing  ~  k 6  nodes. 
Each  connected  component  of  a  group  will  then  correspond  to  a  contiguous  group 
of  faces  in  the  triangulated  Voronoi  diagram.  Using  the  Lipton-Tarjan  planar  sepa¬ 
rator  theorem  [39],  we  can  split  V  in  half  by  removing  0(y/kn)  edges.  This  can  be 
accomplished  in  O(kn)  time.  Continue  to  apply  the  separator  theorem,  at  each  level 
splitting  the  existing  groups  in  half  until  the  point  is  reached  where  groups  contain 
~  k6  nodes.  At  this  poir1  O(p-)  edges  will  have  been  removed  and  we  will  have 
used  0(kn  log(^))  time  running  the  planar  separator  algorithm.  Let  Uj , . . . ,  Vj  be 
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Figure  7-1:  The  triangulated  Voronoi  diagram  V  with  darker  edges  defining  Vi, ,  V/ 

the  resulting  l  =  O(f^)  pieces  of  V.  Then  define  Pt  to  be  the  union  of  the  Sv  such 
that  v  €  Vi. 

Observation  7.4.2  For  any  query  point  q,  the  k-nearest  neighbors  of  q  are  com¬ 
pletely  contained  in  some  P{. 

Proof:  Let  Fq  be  the  face  in  V  which  contains  q  and  vq  be  the  corresponding  node 
in  the  dual  V.  Then  vq  €  V  for  some  i  and  the  set  of  ^-nearest  neighbors  to  q  is 

sVq  c  p,.  □ 

For  all  z,  let  V,  be  the  region  of  the  plane  defined  by  Vi.  See  Figure  7-1.  Notice 
that  collectively  there  are  at  most  O(fr)  connected  sub- regions  among  the  regions 
Vi,...,Vj .  These  connected  sub-regions  are  defined  by  the  O(p-)  edges  of  V  that 
are  removed  during  the  planar  separator  process  and  any  Vi  may  be  the  union  of 
numerous  connected  sub-regions.  Using  Kirkpatrick’s  planar  point  location  algorithm 
[35],  preprocess  this  set  of  edges  to  create  a  data  structure  which  on  input  q,  can 
locate  the  region  V,  containing  q.  Call  this  the  locator  tree  for  V] , . . . , V/ .  This  tree 
has  space  linear  in  the  number  of  edges  needed  to  bound  the  Vi’s.  Furthermore,  since 
Kirkpatrick’s  planar  point  location  data  structure  on  m  regions  can  be  computed 
in  O(miogra)  time,  the  locator  tree  can  be  constructed  in  O(-jplogn)  preprocessing 
time. 

Observation  7.4.3  The  locator  tree  uses  only  O(-jp)  space. 
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Figure  7-2:  Linear  structures  for  the  P,’s  give  a  linear  structure  for  P 

Proof:  Each  boundary  edge  of  a  region  V'  corresponds  to  a  removed  edge  in  the  dual. 
We  have  removed  at  most  O(p-)  edges.  □ 

Once  we  have  located  the  V,  which  contains  q,  we  need  only  search  P,  to  find 
the  A;-nearest  neighbors  of  q.  Since  P,  was  defined  from  ~  k6  sets  of  points  Sv  for 
v  €  Vi,  and  each  Sv  contains  k  points.  ||P,||  <  k7 .  Hence  by  traversing  the  locator 
tree  in  O(logn)  time,  we  can  reduce  the  original  query  problem  to  one  of  finding  the 
^-nearest  neighbors  out  of  a  set  P,  with  only  0(k7)  points.  At  this  point,  one  should 
notice  that  the  sets  Pi,....  Pi  are  not  pairwise  disjoint.  In  fact,  we  have  ~  fz  P,'  s. 
each  of  size  O(k'),  and  a  naive  analysis  would  suggest  that  0(kn )  memory  is  required 
to  store  the  P,’s.  In  order  to  obtain  an  0(n )  upper  bound  on  the  size  of  the  data 
structure  we  are  building,  it  is  necessary  to  give  a  strict  upper  bound  on  the  number 
of  duplications  of  points  which  occur  among  the  P,'s.  Toward  this  end.  we  show  that 
the  total  number  of  points,  counting  duplications,  which  need  to  be  stored  in  the  data 
structure  is  bounded  above  by  ( 1  +  0(£))n: 
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Lemma  7.4.4  £'=1  ||P,||  <  (1  +  0(£))n. 

Proof:  We  begin  by  assigning  each  point  p  £  P  to  the  node  vp  £  V  whose 
corresponding  face  in  V  contains  p.  Now  we  argue  that  the  number  of  duplicates  of 
points  in  P  which  are  contained  in  the  P,’s  is  bounded  by  k  times  the  number  of  edges 
cut  during  the  planar  separator  stage  of  the  construction.  Let  x  €  P,.  Then  x  €  Sv 
for  some  v  £  V{.  So  v  €  ( V)x  which  is  connected  by  Observation  7.4.1.  Then  either 
(V)x  is  completely  contained  in  V,  (and  x  occurs  uniquely  in  P,)  or  else  an  edge  e  of 
(V')x  was  cut  during  the  construction  of  V(  (i.e.,  a  node  of  V{  must  be  incident  to  a 
cut  edge  e  for  which  x  £  Le).  Since  Le  contains  at  most  k  points  for  any  e  and  at 
most  O(p-)  edges  e  have  been  cut,  the  total  number  of  duplicate  points  in  all  of  the 
P.’s  is  0(f).  □ 

7.5  Probabilistic  Techniques 

Each  of  the  sets  of  points  P,  will  now  be  treated  separately.  It  follows  from  Lemma 
7.4.4  that  by  creating  a  linear  size  data  structure  for  ^-nearest  neighbor  queries  on 
each  P,,  we  can  build  a  linear  size  data  structure  for  ^-nearest  neighbor  queries  on 
P  using  the  locator  tree  discussed  above.  See  Figure  7-2  for  a  picture  of  such  a  data 
structure. 

7.5.1  Review  of  Chernoff  Bounds 

Many  of  the  results  in  this  section  use  mathematical  theorems  for  bounding  the  tail 
of  a  binomial  distribution.  These  results  are  referred  to  loosely  as  Chernoff  bounds 
and  can  be  found  in  many  sources  [3,  21,  30,  42]. 

Let  „Yi, . . . ,  Xk  be  independent  0  —  1  random  variables  with  probability  q  of  being 
1.  Then  let  \  =  X,,  and  0  >  1.  We  would  like  to  bound  the  probability  that  \ 

exceeds  its  expected  value  qk  by  more  than  a  factor  0.  The  following  result  is  a  well 
known  Chernoff  bound  and  can  be  easily  derived  using  moment  generating  functions. 
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Lemma  7.5.1  Pr[\  >  0qk]  <  exp  ((0  —  1  —  0  In  0)qk) 

7.5.2  Assigning  the  P7  to  Buckets 

For  each  P,  create  a  set  of  s  =  buckets  B\ , . . . ,  Bf.  We  assign  the  points  of  P, 
to  these  buckets  with  equal  probability.  That  is,  a  point  p  €  P,  is  assigned  to  bucket 
Bj  with  probability  Our  intention  is  to  search  each  bucket  for  the  'N/  log2  k 

nearest  neighbors  of  q ,  and  in  this  manner  capture  the  k-nearest  neighbors  of  q  in 
Pt.  To  do  this,  we  need  to  show  that  for  any  q,  the  k-nearest  neighbors  of  q  are 
evenly  distributed  among  the  buckets.  Consider  the  kth- order  Voronoi  diagram  on  Pt. 
This  diagram  has  m  =  0(k8)  regions  corresponding  to  all  possible  sets  of  k- nearest 
neighbors  to  q.  Call  the  corresponding  sets  of  k  points:  Cj,...,Cm  the  constraints. 
We  say  that  a  constraint  CT  is  satisfied  if  each  bucket  £?/,...,  B ■  contains  less  than 
log2  k  +  0( log  ^  k )  points  from  Cr. 

Lemma  7.5.2  With  probability  >  1  —  L  every  constraint  is  simultaneously  satisfied 
by  a  random  assignment  of  P%  to  the  buckets. 

Proof:  This  is  a  direct  applications  of  the  Chernoff  bound  given  in  Lemma  7.5.1: 
Pr[x  ^  fiqk]  <  exp  (((3  —  1  —  0  In  0)qk)  Choose  a  constraint  CT  =  {pi, . .  .pk}  and  for 
/  =  1  . . .  k  let  Xt  =  1  if  pi  €  Bf  In  this  case  q  =  We  say  that  Cr  is  satisfied  if  less 
that  log2  k  +  7  logi  k  points  from  CT  are  contained  in  any  one  bucket.  We  will  show 
that  7  can  be  chosen  to  satisfy  the  lemma.  Now  let  0  =  1  +  — .  Then  CT  is  satisfied 
if  less  than  0qk  points  from  CT  are  contained  in  any  one  bucket.  Since  there  are  0(ks) 
constraints  and  O(k)  buckets,  from  the  Chernoff  bound  given  above,  we  conclude  that 
the  probability  that  some  constraint  is  not  satisfied  is  less  than  k9Pr[  the  number 
of  points  from  Cr  in  bucket  P-  >  log2  A:  +  7  log^  k ]  <  k9exp((0  —  1  —  0  In  0)qk)  < 
k9exp(  —  7  k  -f  7  Hence,  in  order  to  satisfy  the  lemma,  it  suffices  to  pick  7 

such  that  k9exp(  — 7  -f  7  V^og k |  <  1  When  k  is  sufficiently  large,  7  >  4  is  good 
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enough.  □ 

For  each  set  Pt,  assign  the  points  randomly  to  buckets  and  then  check  all  of  the 
constraints  to  see  that  they  are  satisfied.  Lemma  7.5.2  tells  us  that  with  probability 
>1  —  4  the  assignment  to  Pi  will  check  out  on  the  first  try.  We  can  build  the  kth-order 
Voronoi  rHagram  on  Pi  in  0{k2\\Pi\\  +  jjPt||  log(||P<j|))  time  (l,  36].  By  Lemma  7.5.2  we 
expect  to  check  at  most  2  assignments  in  order  to  successfully  divide  P,  into  buckets. 
Hence  the  expected  time  to  assign  all  of  the  Pi's  to  buckets  is  0(k2n  +  n  log  n).  Here 
we  have  described  a  Las  Vegas  algorithm  to  divide  up  the  regions  P,  into  buckets.  The 
assignments  to  buckets  are  guaranteed  to  satisfy  all  constraints,  but  the  running  time 
of  the  algorithm  is  not  guaranteed.  However,  with  high  probability  this  algorithm 
will  terminate  in  0(k2n)  time.  In  the  next  section  we  show  how  to  assign  the  P,’s  to 
buckets  deterministically. 


7.6  The  Compacted  Data  Structure 

At  this  point,  we  have  a  data  structure  which  divides  P  up  into  sets  Pj, . . . ,  P/  with 
very  little  duplication  of  points.  Given  a  query  point  q ,  we  can  locate  the  proper  set 
Pi  which  contains  q's  ^-nearest  neighbors  in  0( log  n )  time.  Each  Pt  is  further  divided 
up  into  buckets.  We  have  to  retrieve  the  log2  k  +  4  log5  k  nearest  neighbors  to 
q  from  each  bucket  Bj, . . .  ,B*.  The  idea  now  is  to  recursively  build  a  data  structure 
for  retrieving  log2  k  - f  4  log5  k  nearest  neighbors  in  each  bucket.  At  the  bottom  level, 
where  we  are  searching  buckets  for  the  c  (a  fixed  constant  which  does  not  vary  with 
k)  nearest  neighbors  of  q.  we  simply  construct  ct/l-order  Voronoi  diagrams  on  the 
points  in  these  buckets  and  process  these  diagrams  using  Kirkpatrick’s  planar  point 
location  techniques  to  get  linear  space,  0(c)  time  structures  for  finding  the  c-nearest 
neighbors  to  q.  The  c  points  retrieved  from  each  leaf  are  stored  in  a  global  linear  array 
of  possible  neighbors  for  q.  In  what  follows,  we  will  show  that  at  most  0(k)  points  are 
stored  in  this  array.  This  array  of  0(k )  points  is  guaranteed  to  contain  the  Ar-nearest 
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neighbors  of  q.  To  find  the  Ar-nearest  neighbors,  we  will  begin  by  finding  the  kih- 
nearest  point,  x,  to  q  in  this  array  in  O(k)  time  by  using  an  order  statistic  algorithm 
[12].  Then  simply  compare  each  point  to  x  and  keep  those  which  are  closer  or  equal 
in  distance  to  q.  Notice  that  as  the  Ac-nearest  neighbors  algorithm  traverses  the  data 
structure  we  have  constructed  at  the  intermediate  levels  no  points  are  inserted  into 
the  array.  Points  are  only  retrieved  and  written  to  the  array  in  the  leaves,  and  points 
are  not  passed  recursively  up  the  data  structure.  Once  the  entire  structure  has  been 
traversed,  and  all  of  the  relevant  leaves  have  been  investigated,  then  the  ^-nearest 
neighbors  to  the  query  point  are  eliminated  from  among  those  present  in  the  global 
array. 


Below  we  give  a  careful  analysis  of  the  data  structure  to  show  that  it  occupies  0(n) 
space.  We  also  solve  the  necessary  recurrence  to  prove  that  the  structure  reports  a  set 
of  O(k)  points  containing  the  A:-nearest  neighbors  of  a  query  point  q  in  O(logn  +  k) 
time. 

The  data  structure  we  have  constructed  is  a  tree  with  levels  alternating  between 
locator  trees  and  pointers  to  buckets.  Notice  that  at  level  i  of  this  structure  we  are 
breaking  up  a  search  problem  for  k ,  neighbors  into  a  number  of  new  search  problems 

I 

for  ki+i  =  log2  A:,  +  0(\og*  k,)  neighbors.  Hence  at  level  i,  0(log  log  . . .  log  A:)2  is  a 
lower  bound  on  the  range  of  the  neighbor  searches  being  performed.  Now  notice  from 


Observation  7.4.3  that  the  locator  trees  at  level  i  will  occupy  0(J2,  rrM  space,  where 
rrij  are  the  sizes  of  the  buckets  on  level  i  —  1.  By  repeated  applications  of  Lemma 
7.4.4:  mj  <  (1  +  £)  •  •  •  (1  +  ^-L-)n.  Hence,  the  total  amount  of  space  occupied  by 

the  locator  trees  is:  0(p-  +  4-  - — |Qg°^V~ — ♦-•••)  =  0(n )  space.  Since  the 

number  of  bucket  nodes  increases  geometrically  every  other  level,  the  space  used  to 
store  the  bucket  nodes  is  dominated  by  the  space  used  to  store  the  leaves  of  the  data 
structure.  Each  leaf  holds  a  constant  sized  Voronoi  diagram.  Hence  the  total  amount 
of  space  used  to  store  this  structure  is  bounded  by  a  constant  times  the  total  number 
of  points  (including  duplicates)  which  appear  in  the  Voronoi  diagrams  in  the  leaves. 
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By  Lemma  7.4.4  this  number  is  0(((1  4-  ±)(1  4-  1  +  jjsjfejr)-  ■  •)«)  =  0(n). 

The  upper  bound  on  running  time  for  extracting  the  ^-nearest  neighbors  is  equal 
to  the  time  required  to  traverse  the  first  level  locator  tree,  plus  upper  bounds  on  the 
times  required  to  search  recursively  through  the  first  level  buckets: 

T(n,  k)  <  a  log  n  +  - — ^—T(k7,  log2  k  +  A  log?  k) 

log'  k 

T(n,c )  <  ac 


This  can  be  solved  by  substitution  using  the  following  expression  for  T(n,k): 
~  .  6q(1  +  — v4 lalfc  6a(l  4 — 4^)(1  -f-  — 7**— )fc 

6c*k  sj logk '  1 J logk '  v/loglog  k' 

a  log  n  4 - h - - - 4 - - - - - 4.  . .  . 

S  logit  log  log  k  log  log  log  k 


Hence,  we  have  T(n,k)  =  O(logn  4-  fc).In  the  case  that  the  query  point  q  is 
restricted  to  a  set  of  0(n)  points,  then  this  time  can  be  reduced  to  O(fc).  To  see  this, 
let  Ci,...,  Cm  be  the  possible  values  for  q  where  m  =  O(n).  Eliminate  the  locator 
tree  at  the  first  level  of  the  compacted  data  structure  and  replace  it  with  an  array 
A[1 . . .  m]  where  A{j]  contains  a  pointer  to  the  node  on  the  second  level  which  handles 
the  set  Px  containing  the  ^-nearest  neighbors  of  c;.  On  input  Cj,  begin  traversing  the 
compacted  data  structure  at  the  node  indicated  by  the  pointer  .A[/].  In  the  running 
time  analysis  this  replaces  the  alogn  term  with  0(1),  and  the  recurrence  for  T{n.  k) 
now  has  solution  0(k). 

Lastly,  we  must  check  that  at  most  0(k)  points  are  reported.  Notice  that  points 
are  reported  only  in  the  leaves  of  this  structure  and  are  not  passed  up  the  tree  and 
filtered  at  each  level.  In  fact,  if  we  computed  order  statistics  and  filtered  points  at 
each  level,  theu  the  resulting  report  time  would  have  been  fi(logn4-  klog"  k).  Implicit 
in  the  recurrence  given  above  is  the  fact  that  the  superset  of  the  fc-nearest  neighbors 
collectively  reported  by  the  leaves  has  size  at  most:  ((1  4 — 4—  )(1  4 — ,  )(1  4- 

y/logk  y4oglog* 

>7iog  log  log  k )  'H'  =  To  see  this  notice  that  at  the  first  level  we  are  neighbor 
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searching  for  log2  k  +  4  log5  k  points  from  each  of  buckets.  So  at  the  first  level, 
the  number  of  points  which  are  going  to  be  reported  is  bounded  above  by  ( 1  +  —*  )k. 

V  log  k 

A  similar  argument  adds  an  additional  factor  of  (1  + 


\/log  log  k 


)  at  the  second  level. 


etc. 


7.7  Extensions  to  Voronoi  Diagrams  in  Other 
Metric  Spaces 

The  ideas  used  above  to  construct  a  linear-sized  data  structure  for  optimal-time  k- 
nearest  neighbor  searching  are  general  enough  to  apply  to  Voronoi  diagrams  con¬ 
structed  using  a  wide  variety  of  distance  functions.  In  fact,  it  is  clear  from  our  dis¬ 
cussion  that  the  techniques  work  whenever  the  fctA-order  Voronoi  diagram  is  planar 
and  contains  O(k0n )  edges  and  faces  for  some  constant  /3.  The  only  problem  may  be 
that  the  edges  of  a  Voronoi  diagram  in  an  alternative  metric  may  not  be  straight  lines. 
In  that  case,  a  straightforward  application  of  Kirkpatrick’s  planar  point  location  algo¬ 
rithm  will  not  be  sufficient  for  the  construction  of  the  locator  trees.  However,  in  many 
cases  involving  edges  defined  by  curves,  we  may  be  able  to  use  a  planar  point  location 
algorithm  given  by  Edelsbrunner,  Guibas,  and  Stolfi  [28]  and  described  in  Section  7.2. 
Given  that  the  diagram  has  straight  edges  or  that  curved  edges  can  be  handled,  we 
can  always  apply  the  planar  separator  theorem  to  build  a  locator  tree  leading  to  sets 
Pi  containing  0{ka+6)  points.  At  this  point,  we  will  have  0(k2)3+6)  constraints  im¬ 
posed  on  the  buckets  for  Pt.  As  in  the  L2  case,  these  constraints  can  be  satisfied  with 
high  probability.  The  resulting  data  structure  obtained  by  applying  these  techniques 
recursively  satisfies  the  same  recurrences  as  above  with  slightly  different  constants. 
Applying  these  ideas  to  non-euclidean  Voronoi  diagrams  (such  as  those  generated 
by  the  power-distance  metric  or  weighted  metric  with  additive  weights)  is  important 
in  the  applications  presented  in  Chapter  8.  Note  the  the  Voronoi  diagrams  for  these 
metrics  have  k°^  edges  and  the  connectivity  property  described  in  Observation  7.4.1. 
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7.8  Removing  Randomness  from  the  Construc¬ 
tion 

To  make  the  above  construction  deterministic,  we  use  the  well  known  probabilistic 
method  described  in  [42].  For  each  P,,  we  will  proceed  by  first  assigning  the  points  to 
two  buckets  Z?0  and  B\.  Then  each  of  these  buckets  will  be  split  into  two  buckets,  and 
then  into  four  buckets,  etc.,  until  after  log  k  —  2  log  log  k  iterations  the  points  of  P, 
have  been  partitioned  into  buckets.  During  the  initial  iteration,  the  points  of  P, 
are  taken  in  some  order  pi,  p2, .  ■  .  and  assigned  to  B0  or  B The  probabilistic  method 
allows  us  to  ensure  that  each  bucket  ends  up  containing  <  ~  -fi  0(  */k  log  k)  points 
from  each  constraint  Cj.  Central  to  the  methods  in  [42]  for  solving  this  problem  is  the 
computation  of  the  conditional  probabilities:  Pr[  >  |  +  4 C k  log  k  points  from  Cj  are 
assigned  to  B0  |  px  6  B0,p2  €  Bx, . . .  ,pm  £  B0  ].  Call  this  conditional  probability  Pr[ 
Cj  |  Pi  €  B0,  p2  €  jBi  , . .  ■ ,  Pm  €  Bo  ]•  Such  a  conditional  probability  can  be  computed 
exactly  as  a  binomial  series  with  <  k  terms.  We  will  come  back  to  the  problem  of 
computing  the  Pr[  Cj  |  px  €  B0,p2  €  Si, . . .  ,pm  €  S0  ]’s  shortly.  First,  however,  we 
detail  the  probabilistic  method  used  to  assign  the  points. 

Suppose  that  points  pi, . .  .  pm  have  been  assigned  to  buckets  So  and  B\  in  such  a 
manner  that  there  still  exists  an  assignment  of  the  remaining  points  which  does  not 
violate  any  of  the  constraints.  We  would  now  like  to  assign  pm+J  to  one  of  the  buckets 
in  a  manner  which  preserves  the  property  that  an  assignment  of  the  remaining  points 
exists  which  does  not  violate  any  of  the  constraints.  So  for  each  of  the  constraints. 
Cj,  effective  in  set  Pt,  we  compute  the  conditional  probability  for  B  =  B0,  B\ :  Pr[  C: 

I  Pi  €  B0,p2  €  pm  €  B0 ,  pm+i  €  B}. 

By  hypothesis,  for  either  B  =  B0  or  B  =  B\  (or  perhaps  both)  the  sum  of  these 
conditional  probabilities  will  be  less  than  1:  Pr[  !  Pi  €  B0,p2  €  Bx,....pm  € 

Bo  Pm+i  €  B\  <  1.  To  see  this,  notice  that  we  have  assumed  as  our  induction 
hypothesis  that  Ylj  Pr(  Q  |  Pi  €  B0,P2  €  Bx, pm  €  f?0]  <  1-  But  Pr[  C;  | 
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P\  G  B0,  P2  G  B\, . . . .  pm  G  B0]  —  ^Pr(  Cj  |  pi  G  Bo.  P2  G  Bi , . .  . ,  pm  G  Bo  pm+ 1  €  B0] 

4-  iPr[  Cj  |  P!  €  B0,p2  G  Bx . pm  G  B0  pm+ 1  G  Bj].  Hence,  |  Pr[  Cj  |  pi  G 

Bq,P2  G  B\, - pm  G  Bq  pm+i  G  B0]  +  jUj  Pr[  Cj  |  p\  G  Bo,p2  G  Bj . . . . . pm  G  Bo 

pm+1  G  Bi]  <1.  It  follows  that  one  of  the  two  sums  on  the  left  hand  side  of  this 
inequality  must  be  less  than  1.  In  the  base  case,  where  none  of  the  points  have  been 
assigned  to  any  buckets  yet,  Lemma  7.5.2  assures  us  that  the  initial  sum  will  be  less 
than  1.  We  then  assign  pm+i  to  a  bucket  which  achieves  a  favorable  set  of  conditional 
probabilities.  In  the  course  of  this  process,  we  have  guaranteed  that  each  bucket  gets 
split  roughly  in  half  for  log  k  —  log  log  k  levels.  By  construction,  the  maximum  number 
of  points  from  any  constraint  in  a  given  bucket  at  stage  i  obeys  the  following 
recurrence:  m,  <  "T"1  +  0( y/mi~[JogT)  where  mo  =  k.  It  follows  by  substitution 
that  m,  <  £  +  0(^f  £  log  k)  for  all  i  such  that  -E  =  ff(logfc).  Hence,  at  the  stage 
numbered  log(i^pr^)i  we  have  produced  roughly  equal  sized  buckets  Pi,...  P  ^  such 
that  each  bucket  contains  at  most  log2  k  4-  0( log  ky/IogTc)  points  from  any  C:. 


We  now  turn  our  attention  to  computing  the  conditional  probabilities.  Let  u  be 
the  number  of  points  from  C j  which  have  not  been  assigned  to  buckets  yet.  Let 
v  =  +  4\/k  log  k)  -  ^(points  from  Cj  that  have  been  assigned  to  B0).  Define  c(u.  u) 


—  _L 

2U  t_i=v 


Then  P1"]  Cj  |  pi  G  B0,p2  G  B1? - pm  G  B0  ]  =  c(u,u). 


To  compute  these  conditional  probabilities  efficiently,  it  will  suffice  to  precompute 
(  U  \  for  u  <  k  and  /  <  k.  Since  f  “  +  1  )  =  -S±*t  (  “  )  and  (  “  \  = 


,  this  precomputation  can  be  achieved  iteratively  in  0(k 2)  time.  Iteratively 


precompute  all  of  the  c(u.  v)  for  all  u  <  k.  -  <  t>  <  k.  We  can  then  simply  look  up  the 
conditional  probabilities,  and  each  Pr[  Cj  |  pi  G  B0,p2  G  Bi,...,pm  G  B0  ]  requires 
0(1)  time  to  lookup.  The  number  of  constraints  active  in  B,  is  bounded  above  by  the 
number  of  faces  in  the  kth -order  Voronoi  diagram  on  B,:  0(Aj|BIj|).  Hence,  there  are 
0(&||B,||)  active  constraints.  Cj,  in  B,.  After  precomputing  all  of  the  c(u.v),  Pr{ 


so 
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Cj  |  pi  €  B0-P2  €  B 1, . . .  .pm  6  Bo)  requires  0(A:j|P,||)  to  compute.  This  computation 
is  performed  0(\ogk )  times  to  assign  all  points  to  buckets.  Hence,  the  total  time 
required  to  assign  points  to  buckets  is  0(kiogk\\P,\\  -f  k2)  =  0(k  log  A'|jP,  j|).  Hence, 
the  entire  assignment  process  for  the  first  level  of  the  structure  takes  0{k  log  kn )  time. 
At  subsequent  levels,  we  are  building  structures  for  the  fc, -nearest  neighbor  problem, 
where  kt  =  log  kt  +  0( log1  kx)  Hence,  the  bound  0{k{  log  fc,n)  for  assigning  points  to 
buckets  at  each  level  decreases  more  than  geometrically  fast. 


7.9  Total  Preprocessing  Time 

The  total  computation  time  for  building  the  compressed  Voronoi  diagram  is  given  by 
the  sum  of  the  times  for  the  following  operations.  The  time  given  in  each  case  is  the 
sum  of  the  associated  operations  performed  at  every  level: 

1.  compute  all  intermediary  Voronoi  diagrams:  0(k2n  +  nlogn)  [1,  36]. 

2.  compute  planar  separators  to  build  Px' s:  (kn  log-^-)  [39]. 

3.  build  all  of  the  locator  trees:  O(n)ogn)  [35]. 

4.  A  deterministic  assignment  of  the  Pi's  to  buckets:  O(k2nlogk)  [Section  7.S]. 

5.  A  randomized  assignment  of  the  P,'s  to  buckets:  0(k2n  4-  n  log  n)  expected  time 
’Section  7.5]. 

At  the  ith  level  of  the  construction,  we  are  dealing  with  &,-nearest  neighbor 
problems  where  k\  =  k  and  ki+i  =  log  &  +  log5  kt.  (For  convenience  we  desig¬ 
nate  k0  =  n.)  Since  the  kth  order  Voronoi  diagram  can  be  computed  in  0(k2n  + 
n  log  n  j  time,  at  the  ith  level  computing  the  Voronoi  diagrams  for  all  the  Pj s  costs 
0(kj]\P ,||  +  |j/j||  log  ||/>||)  =  Y.P,  e(^l|e,||  +||Cj||  log  )  =  0(/:?n  +  nlogfr,-i). 

Summing  over  all  levels,  we  see  that  the  total  time  spent  computing  Voronoi  diagrams 
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is  0(k2n  +  n  logn).  Likewise  as  we  showed  in  Section  7.4,  at  level  i  the  planar  separa¬ 
tors  computed  to  break  up  Pj  require  0(fc,||Pj||  log(^^))  time  to  compute.  Summing 
over  all  i  gives  the  bound  stated  above.  Finally,  we  have  shown  above  that  the  locator 
trees  at  level  i  take  time  O(-^logn)  time  to  construct.  Summing  over  all  i  gives  a 
total  preprocessing  time  of  0(n  logn)  for  building  locator  trees. 

The  dominating  term  in  the  sum  of  these  time  bounds  varies  with  k,  but  it  is 
clear  that  the  construction  can  always  be  accomplished  in  O(k2n\ogk  +  knlogn) 
deterministic  time  or  0(k2n  +  kn\ogn)  randomized  time. 


7.10  A  Monte  Carlo  Construction 

Using  a  monte  carlo  construction,  we  can  build  a  compacted  data  structure  with 
0(n  log2  n  log  log  n)  preprocessing  time.  For  k  <  log  n,  the  deterministic  construction 
has  preprocessing  time  0(nlog2n).  For  k  >  logn,  the  preprocessing  time  for  the 
above  deterministic  construction  is  dominated  by  the  computation  of  Voronoi  dia¬ 
grams  and  the  assignment  of  points  to  buckets  at  the  one  or  two  top  levels.  However, 
when  k  >  logn,  we  can  discard  the  initial  computation  of  Voronoi  diagrams,  and 
assign  points  directly  to  buckets  in  a  randomized  manner.  In  what  follows  we  show 
that  with  high  probability  this  approach  builds  a  correct  data  structure  for  solving 
the  fc-nearest  neighbors  problem  in  0(n  log2  n  log  log  n)  preprocessing  time.  Further¬ 
more,  each  time  we  query  this  data  structure  we  get  a  witness  which  tells  us  whether 
or  not  the  answer  provided  to  our  query  is  correct. 

When  k  >  logn,  skip  the  initial  planar  separator  stage,  and  at  the  first  level 
assign  the  points  directly  to  buckets  using  0[n)  time.  As  we  demonstrated  in 
the  previous  analysis  of  probabilistic  bucket  assignment,  with  high  probability,  in  or¬ 
der  to  capture  the  fc-nearest  neighbors  of  q  it  suffices  to  report  the  0{\ogn)  nearest 
neighbors  of  q  in  each  bucket.  This  can  be  accomplished  using  the  above  determinis¬ 
tically  constructed  data  structures  requiring  0(n  log2  n  log  log  n)  total  preprocessing 


82 


CHAPTER  7.  VORONOI  DIAGRAM  COMPACTIO.X 


time.  Hence  the  total  preprocessing  time  to  build  this  Monte  Carlo  structure  is 
0(n  log2  n  log  log  n).  If  we  use  a  randomized  assignment  of  points  with  checking  (as 
presented  in  Section  7.5)  at  the  first  level,  this  time  bound  then  reduces  to  0(n  log2  n). 

Given  a  query  q  to  this  new  data  structure,  keep  track  of  the  O(logn)  points 
reported  from  each  of  the  first  level  buckets.  We  are  certain  that  each  bucket  reported 

correctly.  Now  let  p j, - pk  be  the  final  report  from  this  new  structure.  If  our  top 

level  assignment  of  points  to  buckets  is  correct,  there  should  be  at  least  one  point 
reported  from  each  bucket  that  does  not  make  the  final  list  pl,...,p/e.  These 
discarded  points  are  our  witnesses  that  the  final  list  is  correct.  If  all  of  the  points 
reported  from  some  bucket  are  contained  in  the  final  list,  then  that  bucket  contains 
more  than  log  k  +  4  logs  k  of  the  nearest  neighbors  to  q  and  the  top  level  assignment 
of  points  to  buckets  violates  at  least  one  constraint. 


Chapter  8 

Applications  of  the  Compaction 
Technique 


8.1  k- Nearest  Neighbor  and  Circular  Range 

Search 

Given  a  set  of  n  points  in  the  E'ididean  plane  E 2 .  and  a  fixed  k.  the  technique  given 
in  Chapter  7  provides  an  0(n)  space  data  structure  for  computing  the  ^-nearest 
neighbors  of  a  query  point  q.  However,  if  k  is  not  fixed,  but  rather  given  as  a  part 
of  the  query,  then  we  simply  construct  log  n  —  log  log  n  linear  size  data  structures 

D] . Du  where  the  D,  can  be  used  to  obtain  the  2’ log  n- nearest  neighbors  of  q. 

These  0( log  n )  structures  occupy  a  total  of  0(n  log  n )  space.  Now,  given  a  pair  ( q.  k). 
we  simply  compute  the  smallest  j  such  that  2J  >  L^j  and  use  Dj  to  compute  the 
2J  log  n- nearest  neighbors  of  q  in  O(logn  4-  2;  log  n )  =  O(logn  ~f  k)  time.  From  this 
set  of  0(k  log  n )  neighbors,  filter  out  the  ^--nearest  in  0(k  +  logn)  additional  time 
by  finding  the  ktk  order  statistic  i  [12]  and  keeping  the  points  which  are  as  close  to 
q  as  x. 


To  solve  the  circular  range  search  problem  we  also  use  the  logn  —  log  log  n  data 
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structures  described  above.  When  given  a  query  disk  with  center  c  n nd  rad:us  r,  use 
q  to  query  D\ ,  D 2, . . .  until  some  structure  D:  provides  output  which  contains  a  point 
x  farther  than  distance  r  from  q.  This  is  a  straightforward  application  of  the  filtering 
search  technique  introduced  by  Chazelle  [13]. 

The  central  idea  in  filtering  search  is  that  the  search  and  report  parts  of  a  query- 
retrieval  algorithm  should  be  made  dependent  upon  each  other.  In  this  manner,  the 
more  points  the  algorithm  is  going  to  report,  the  longer  it  is  allowed  for  searching. 
We  can  see  that  considering  Di,D2,-..  in  order  is  an  effective  means  of  balancing 
searching  and  reporting  in  the  circular  range  search  problem.  If,  after  searching  Dt. 
we  find  that  all  points  reported  are  still  closer  than  distance  r  from  point  q,  we  then 
know  that  the  report  part  of  the  algorithm  will  run  for  at  least  ff(2')  time  (to  report 
H(2‘)  points)  and  hence  we  have  enough  time  to  search  Dl+\  without  degrading  the 
asymptotic  running  time  of  the  algorithm. 

Now,  let  j  be  the  smallest  index  such  that  D}  outputs  a  point  1  outside  the  query 
disk.  Then  the  time  to  execute  these  j  stages  is  simply:  O(Y^=0(2'  +  l)logn)  = 
0(2-'  logn  +  logn)  =  0(k  +  logn)  since  k  =  0(2J  log n).  Finally,  search  through  the 
2Jlogn  points  given  by  and  output  those  which  are  contained  inside  the  query 
disk. 


8.2  Half  Space  Range  Search  and  Power  Dia¬ 
grams 

Recall  that  for  the  half-space  range  search  problem  in  3  dimensions,  we  are  given  a 
set  P  of  n  points  in  £3  and  a  query  half-space.  We  are  asked  to  report  all  points  of 
P  that  lie  inside  (or  on  the  boundary  of)  this  half-space.  In  this  section,  we  show 
how  this  problem  can  be  transformed  into  a  nearest  neighbor  search  problem  in  the 
power-distance  metric.  To  begin,  we  show  how  to  transform  the  range  search  problem 
into  a  ray-stabbing  problem  using  geometric  duality  [19,  29]. 
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Geometric  Duality 

In  the  half-space  range  searching  problem,  we  are  given  a  set  of  points  in  3-space 
(a,,  bi,  c,)  for  i  =  1  . . .  n.  Then  for  each  query  half-space  ax  +  0y  +  72  >  1  we  are 
asked  to  return  those  points  which  are  contained  in  the  half-space  or  on  its  boundary. 
Geometric  duality  defines  that  each  point  (a,-,  c,j  has  a  dual  half-space  given  by 
the  equation  a,x  +  bty  +  CiZ  >  1.  Similarly,  the  query  half-space  has  a  dual  point 
(a,  0, 7).  Then  we  have  the  following  duality  result: 

Observation  8.2.1  A  point  (a,,  6,,  c,)  lies  inside  the  half-space  ax  +  3y  +  72  >  1  iff 
the  dual  point  (a.3,~f)  lies  inside  the  dual  half-space  a,x  +  b,y  +  c,2  >  1. 

Proof:  Both  sides  of  the  “iff”  are  equivalent  to:  oa,  -f  /?6,  4-  7 Cj  >1.  □ 

Such  simple  duality  results  have  powerful  consequence  in  computational  geometry. 
[19]  For  our  purposes.  Observation  8.2.1  indicates  that  the  transformation  from  a  half¬ 
space  range  searching  problem  to  a  ray  stabbing  problem  can  be  accomplished  in  0(n) 
time  and  space.  After  the  transformation,  we  are  given  a  set  of  half-spaces  defined  by 
a,x  +  b,y  +  CiZ  >  1  for  i  =  1 ...  n  and  a  query  point  (a,  0,  7).  Our  algorithm  should 
report  those  half-spaces  which  contain  this  query-point. 

Without  loss  of  generality,  rotate  the  coordinate  axes  so  that  none  of  the  planes 
defining  our  half-spaces  is  parallel  to  the  x-axis.  Each  half-space  now  intersects  the 
x-axis.  Define  those  half-spaces  which  contain  a  segment  of  the  z-axis  heading  for 
+00  to  be  upwardly  oriented.  Define  the  others  to  be  downwardly  oriented.  We  will 
split  these  groups  into  two  sets  and  treat  them  separately.  This  will  cause  us  to 
run  the  algorithm  twice,  once  on  the  upwardly  oriented  half-spaces  and  once  on  the 
downwardly  oriented  half  spaces.  Within  a  constant  factor,  however,  this  will  not 
affect  our  running  time.  So.  without  loss  of  generality,  we  will  assume  that  all  of 
the  half-spaces  are  upwardly  oriented.  Hence,  if  we  consider  a  ray  which  originates 
at  point  {a.  3. 7)  and  extends  downward  to  2  =  —00  running  parallel  to  the  z-axis. 
the  point  (a.  3. 7)  is  contained  in  exactly  those  half-spaces  whose  boundary  planes 
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the  ray  intersects.  In  this  manner,  we  have  transformed  the  half-space  range  query 
problem  into  a  vertical  ray-stabbing  problem. 

In  order  to  describe  how  to  solve  the  ray-stabbing  problem  using  &(/!-order  power 
diagrams  we  need  to  understand  more  about  the  relationship  between  power  diagrams 
and  planes  in  euclidean  3-space. 

Power  Diagrams  and  Planes  in  3-Space 

Aurenhammer  and  Edelsbrunner  [6,  5,  27]  provide  a  correspondence  between  the 
arrangements  of  planes  in  3  dimensions  and  power  diagrams  in  2  dimensions.  This 
correspondence  hinges  on  the  notion  of  a  k- set.  Given  a  set  of  n  planes  in  E 3,  the 
k-set  is  the  locus  of  all  points  in  E 3  that  satisfy  the  following  properties: 

1.  Each  point  in  the  k- set  belongs  to  one  of  the  n  planes. 

2.  For  every  point  in  the  k-set  there  are  exactly  n  —  k  planes  passing  above  it  in 
the  vertical  dimension  (i.e.,  the  positive  2-direction). 

Two  points  are  said  to  be  in  the  same  face  of  a  k-set  if  they  lie  in  the  same  plane  and 
have  exactly  the  same  planes  passing  above  them. 

In  the  above  definition  of  a  k- set.  we  say  that  a  plane  lies  above  a  point  (a.  6,  c), 
if  for  some  positive  A,  the  point  (a,  6,  c  4-  A)  belongs  to  the  plane.  Notice  that  the 
boundary  of  a  k- set  is  formed  by  the  lines  defining  the  intersection  of  pairs  of  planes. 
Let  A  be  a  face  in  a  Ar-set.  Then  we  define  the  projection  of  A  to  be  the  image  of 
A  projected  onto  a  plane  at  2  =  —00:  {(x,y):  ( x,y,z )  €  A  for  some  2}.  Now  fix- 
some  k0.  It  is  not  hard  to  see  that  the  projection  of  all  the  faces  of  the  k0- set  gives  a 
division  of  the  plane  into  convex  sub-regions. 

We  now  wish  to  establish  a  correspondence  between  the  weighted  points  which 
define  a  power  diagram,  and  planes  in  3-space.  For  proofs  and  more  detailed  ex¬ 
planations  of  the  results  appearing  beiow\  we  refer  the  reader  to  the  work  of  Au¬ 
renhammer  and  Edelsbrunner  [6,  5,  26.  27] .  Recall  that  each  point  p,  in  the  set  of 
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points  defining  a  power  diagram  has  an  associated  weight:  tt’(p,).  Define  a  plane 
7 r(p,):  2  =  2(x,y)Tpx  —  pf p,  +  w(px)  where  [x,y)Tpx  is  the  dot  product  of  (x,y)  and 
(x,-, yx)  and  p,  =  (x,-,y,).  Furthermore,  for  any  two  points  pxxp3  we  define  chor(p,, p3) 
to  be  the  locus  of  points  equidistant  from  p,  and  p3  in  the  power  distance  metric 
with  weights  w(px)  and  u>(pj).  It  is  not  hard  to  see  that  chor(p,,  p})  is  a  straight  line. 
Furthermore,  we  have  the  following  observation: 

Observation  8.2.2  The  projection  of  ^{pi)f]^{Pj)  onto  the  plane  is  chor(px,pJ). 

In  fact,  a  much  stronger  result  can  be  proven.  Consider  the  k- sets  formed  by  the 
set  of  planes  T  =  7r(pi),  7r(p2), . . .  ~{pn)-  Then  for  a  fixed  k  we  have: 

Theorem  8.2.3  The  projection  of  the  faces  of  the  k-set  of  T  onto  the  plane  at 
z  =  — oo  is  the  kth -order  power  diagram  on  the  points  px, . .  .  ,p„  with  weights  w(px). 
Furthermore,  if  v(rx)  is  the  plane  containing  one  of  the  faces  A  of  the  k-set,  and 
7 r(r2), . . . ,  7r (rjt)  are  the  planes  below  A,  then  the  k-nearest  neighbors  to  the  points  in 
the  projection  of  A  are  rx, . . . ,  r*. 

Theorem  8.2.3  can  be  used  to  give  an  upper  bound  on  the  size  of  the  Arf/l-order 
power  diagram.  Since  the  maximum  number  of  faces  in  a  k- set  in  3  dimensions  has 
been  shown  to  be  0{k2n )  (see  Clarkson  and  Shor  [24]),  the  number  of  edges  in  a 
kth-ordeT  power  diagram  is  0{k2n).  Consequently,  if  k  is  known  ahead  of  time,  then 
our  compaction  techniques  can  be  used  to  obtain  an  0(n )  data  structure  which  can 
be  used  to  compute  the  /^-nearest  neighbors  in  the  power  metric  in  Offc  +  logn)  time. 
On  the  other  hand,  if  k  is  not  fixed,  but  given  as  a  part  of  the  query,  then  we  can 
construct  logn  —  log  log  n  linear  data  structures  as  in  the  previous  section  and  obtain 
the  fc-nearest  neighbors  of  q  in  O(logn  +  k)  time.  Again,  this  structure  requires 
0(n  log  n)  space. 

Returning  to  the  vertical  ray-stabbing  problem,  px _ _ pn  and  w(px ),....  w{pn) 

are  the  points  and  weights  corresponding  to  the  power  diagram  generated  by  the 
projection  of  f.  We  are  given  a  query  point  (a.  3 , 7),  and  are  asked  to  determine  the 
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planes  that  this  point  lies  above.  Let  Dx,  D2, . . . ,  Aogn-iogiogn  be  the  compacted  data 
structures  described  above.  First,  we  query  D\  with  the  point  ( a,/3 )  to  determine  the 
logn-nearest  points  77, . . . ,  riogn  in  the  power-distance  metric.  This  takes  O(logn) 
time.  Then  tr(ri), . . . ,  tr(riogn)  are  the  lowest  logn  planes  which  might  appear  below 
(a,/3.  We  can  then  check  each  of  these  planes  tt (r\)  in  constant  time  to  see  if 
(a,  /?.  * )  appears  above  7r(r,-).  At  this  point,  we  have  still  used  0(log  n)  time.  If  all  of 

7r(ri) _ ,tr(riogn)  appear  below  (a,/3,7),  then  continue  on  to  check  D2  and  retrieve 

the  2 logn  nearest  neighbors  to  (a,/3).  Proceed  in  this  manner,  using  the  filtering 
search  techniques  presented  in  the  previous  section.  At  some  Dt,  we  will  find  the  first 
plane  r(r_,)  which  lies  above  (a,/3, 7).  At  this  point,  stop  and  report  all  of  the  planes 
from  the  set  77(77), ...,  7r(riogn2.)  which  lie  below  (q,/9,  7).  Let  us  say  that  k  planes 
are  reported.  Then  the  set  of  planes  retrieved  from  Dx  has  size  at  most  2k.  Since  at 
a  minimum  Dx  gives  us  logn  planes  to  search  through,  it  is  clear  that  this  algorithm 
runs  in  O(\og  n  +  k)  time. 


8.3  A:-Nearest  Neighbors  in  the  Weighted  Metric 

Finally,  ^-nearest  neighbor  searching  in  the  weighted  metric  (with  additive  weights) 
can  also  be  performed  using  data  structures  built  from  the  compaction  techniques 
given  in  this  paper.  As  described  in  Section  7.1.2,  the  kth-o rder  Voronoi  diagram  for 
this  metric  has  only  O(kn)  edges.  However,  since  these  edges  are  defined  by  second 
degree  curves  we  cannot  use  Kirkpatrick’s  planar  point  location  algorithm  to  build 
the  locator  trees  in  the  compact  data  structure.  In  this  case,  we  use  the  planar  point 
location  data  structure  given  by  Edelsbrunner,  Guibas,  and  Stolfi  [28]  and  described 
in  Section  7.2.  With  this  change,  the  resulting  compact  Voronoi  diagrams  in  the 
additive- weight  metric  give  the  same  optimal  time  and  space  complexity  performance 
for  finding  fc-nearest  neighbors  as  the  euclidean  structures.  When  k  is  fixed,  the  data 
structure  occupies  0(n )  space  and  responds  to  queries  in  0(\ogn  -|-  k )  time. 
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8.3.1  Further  Extensions  of  the  Compaction  Technique 

The  technique  given  in  Section  7  for  the  ^-nearest  neighbors  problem  may  have  ap¬ 
plications  to  other  retrieval  problems.  For  example,  this  technique  can  be  modified 
to  obtain  an  0(n)  space  data  structure  that  can  be  used  to  solve  a  retrieval  problem 
in  Q(k  +  logrc)  time  for  any  class  of  objects  with  the  following  properties: 

1.  Each  object  in  the  class  is  a  bounded  or  unbounded  region  of  the  plane  and  its 
boundary  is  some  Jordan  curve. 

2.  If  we  consider  any  n  objects  in  this  class  then  the  number  of  disjoint  k-regions  is 
O{k0n)  (for  some  constant  0  >  0).  We  define  a  ^-region  as  a  connected  region 
of  the  plane  such  that  any  point  contained  in  that  region  belongs  to  exactly  k 
objects.  Furthermore,  each  k- region  must  be  adjacent  to  at  most  0{ka)  other 
regions  for  some  constant  a  >  0. 

3.  The  partitioning  of  the  plane  into  ^-regions  by  the  set  of  Jordan  curves  J 
defining  the  n  objects  must  satisfy  certain  point  location  properties:  given  any 
sub-partitioning  of  the  plane  by  some  0(m)  sized  subset  of  J,  there  exists  a  data 
structure  of  O(m)  size  which  can  perform  planar  point  location  in  0( log  m)  time. 

It  is  not  hard  to  see  that  most  geometric  objects  (circles,  ellipses,  polygonal  curves 
with  at  most  a  constant  number  of  edges)  satisfy  these  properties.  Consequently,  for 
most  of  these  problems,  given  a  fixed  value  of  k,  we  can  obtain  an  0(n)  data  structure 
that  reports  the  Q(k)  objects  containing  a  query  point  q  in  0(logn  -f  k)  time.  (If 
more  or  less  than  O(fc)  objects  contain  q  then  the  structure  would  report  the  empty 
set  in  (3(logn)  time.)  However,  for  all  such  cases,  these  results  are  not  new  since  the 
techniques  given  in  [16,  17]  can  be  used  to  obtain  simpler  0(n)  space  data  structures. 
In  a  similar  vein,  we  can  use  this  technique  together  with  filtering  search  to  solve  the 
3-dimensional  dominance  reporting  problem  in  0(n  log  n)  space  and  0(logn  +  £)  time, 
but  again  a  simpler  data  structure  that  achieves  the  same  space  and  time  bounds  has 
been  given  by  Gabow  et  al.  [31]. 
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